Prediction of Maximum Scour around Circular Bridge Piers using Semi-Empirical and Machine Learning Models

1. Introduction

Scour, a natural occurrence, results from the erosive force of flowing streams on erodible beds, risking the integrity of the bridge structure [1,2]. The study of bridge scour has been divided into two major fields scientific inquiry seeks to understand the scour process academically, while engineering endeavors aim to offer practical guidance for implementation such as scour depth (d_s) estimating (flume tests, empirical formula, data-driven soft computing technique) monitoring and countermeasure [3,4].

Many scour prediction formulas have been developed to predict equilibrium scour depth (d_se) using flume tests [5,6,7,8,9], temporal scour variations [10,11,12,13], and maximum scour depth (d_sm) incorporating the time factor [14,15]. Several distinguished researchers have contributed significantly to the design method for time-dependent local scour as well [16,17,18]. After investigations by several researchers, an understanding of d_sm around circular piers has been achieved mostly for sand beds [19,20] and very few for gravel beds [21] in clearwater conditions. The data-driven approach using different ensemble machine learning models to address circular bridge piers [22,23,24,25,26,27] has been studied recently. Empirical formulas for live-bed and clearwater scour were both developed by several researchers [28,29,30,31].

Empirical formulas: Major findings and limitations of some literatures are addressed below. Lacey [32] and Inglis [33] developed regime formulas using data from irrigation canals in India. Kothyari et al. [34] commented on the Lacey–Inglis method specifically designed for non-cohesive sands with sediment mean size (D₅₀) of 0.15-0.43 mm and gradation (σ) of 1.4-1.8, used for designing railways and road bridges in India [29]. Vijayasree and Eldho [35] found that the IRC formula underestimates d_s by overlooking pier geometry. They modified IRC, with an empirical constant from experiments on various pier shapes. Gao et al. [36] developed a formula from field data, applied in Chinese infrastructure for over 20 years. This study refined the formula using hydraulic models and field data to closely align calculated results with real-world observations for scour scenarios. Ansari and Qadar [37] tested three formulas [38,39,40] for predicting d_s, analyzing 100+ field observations from 12 sources and regions. Melville [41] presented a method to predict d_s at piers using laboratory data. However, variations in formula structures and sensitivity of parameters can lead to significantly different predictions [42]. Federal Highway Administration guidelines recommended estimating design d_sm [7]. It was developed from lab data for circular piers by Chabert and Engeldinger [43] and Colorado State University data, and widely acknowledged as a leading formula for predicting d_sm. Benedict and Caldwell [44] reviewed literatures to find pier scour, compiling 569 lab and 1,858 field measurements. Sheppard et al. [45] examined 17 formulas with 441 lab and 791 field data points. They ranked the formulas statistically and compared under-prediction errors to total errors. Manes and Brocchini [46] developed a scale-independent formula for d_se at a circular pier in a granular bed. Vonkeman and Basson [47] evaluated 30 empirical formulas with their 48 experiments and revealed different results, especially for non-circular piers. Coscarella et al. [48] presented new hypotheses on momentum transport and bed shear stress, questioning conventional thinking after Manes and Brocchini [46], and developed a formula for predicting d_sm. Franzetti et al. [14] developed another d_sm formula integrating metadata from literature. In addressing temporal d_s measurements, Nandi and Das [15] compared 11 literature formulas and developed one formula to predict d_st. Equilibrium methods for predicting d_s face well-documented limitations, with studies noting inconsistencies in d_se and t_e definitions, which lead to ambiguity and inaccuracies in predictions [49,50]. This lack of a universal criterion complicates experiments, field studies, and data analysis, contributing to database inconsistency and uncertainty [51,52]. The evaluation of 74 empirical prediction formulas (Section S1) illustrates the varied strengths, limitations, and drawbacks of these methods (Section S2). Many existing formulas for predicting d_s used limited parameter sets based on specific experimental conditions, which may introduce some degree of inaccuracy [53]. While some formulas tend to show good agreement within the datasets they were originally derived from [54,55]. Some formulas may slightly overestimate or underestimate d_sm compared to observed data, highlighting areas where predictive accuracy could be further refined [56,57]. The detailed dimensional information comparing the thorough literature review of experimental and numerical studies around isolated (Table S1) and complex piers (Table S2) is given respectively in the supplementary information. Some approaches excel under specific conditions but may be less accurate across varied datasets, suggesting that a broader parameter inclusion might enhance their applicability.

Data-driven stand-alone and hybrid ML models: The studies have used GEP techniques to predict d_s around bridge piers, significantly improving prediction accuracy compared to traditional empirical models [58]. Rathod et al. [59] demonstrated that support vector machine (SVM) outperformed other artificial intelligence algorithms like ANN, adaptive neuro-fuzzy inference system (ANFIS), and gene expression programming (GEP) in predicting d_s. Eini et al. [60] employed hybrid models combining particle swarm optimization and XGBoost, finding that RPSO-XGBoost provided better results for d_s estimation. Kumar et al. [25] utilized ensemble methods like bagging regressor and adaboost regressor, showing enhanced performance over individual models like SVR. Shalini et al. [61] highlighted the success of the ANFIS model, outperforming other techniques, while Ahmadianfar et al. [62] improved d_s predictions around pile groups using an ensemble approach. Pandey et al. [63] found that the CatBoost model outperforms extra tree regression and K-nearest neighbor in predicting d_s around abutments. Choudhary et al. [64] developed ANFIS and GEP models to accurately predict d_s around bridge piers under clear-water and live-bed conditions. Using key hydraulic and sediment parameters, the ANFIS model proved highly accurate (CD ≈ 0.95) across most cases, though it showed limitations when the pier was relatively narrow compared to flow depth in live-bed scenarios. Recent work by Baranwal and Das [65] and Kumar et al. [66] continued to explore the efficacy of models like M5Tree, XGBoost, and random forest, demonstrating superior predictive accuracy across different datasets. Niknam et al. [67] showed that GEP with a three-gene configuration provided the most accurate predictions for d_s around oblong piers, outperforming SVM and ANN models. A detailed table summarizing the studies that employed GEP and ANN as at least one of their ML methods for isolated pier d_s modeling was included as Table S3 in the supplementary information. GEP and ANN are particularly effective for predicting d_s because they are relatively easy to use with available software frameworks. These models were chosen for simplicity and strong predictive abilities, making them practical without knowing advanced programming skills.

The present research, therefore, aims to formulate robust data-driven formulas using ANN, GEP, semi-empirical approach to analyze. The present study uses 40 experimental datasets, including 768 data points obtained from 67 years of literatures. The number of datasets considered is the highest compared to previous studies. Eight functional parameters are used for inter-comparison of 74 literature formulas (highest so far), including flow intensity(V/V_c), flow shallowness (H/B), sediment gradation (σ), sediment coarseness (B/D₅₀), dimensionless time (Vt/BΔ^0.5), constriction ratio (B/W), Reynolds number (R) and Froude number (F). To capture the relationship between these parameters, different nonlinear functions were used, such as two parameters sigmoid function, bragg function, power function, exponential function, etc. The Gauss-Newton method and Levenberg-Marquardt algorithm are used to optimize the objective function. New experiments are conducted to assess the prediction performance of formulas.

This study offers the following key contributions (1) utilizes a comprehensive dataset of 768 data points collected from various experimental studies, ensuring robust model training/calibration and validation, (2) improving the representation of dimensionless parameter effects in a semi-empirical framework through novel functional formulations, (3) including important factors like B/W and F, which most earlier models overlooked, to better predict d_sm under different flow conditions, (4) applying ANN and GEP techniques to an extensive, d_sm dataset, to improve the prediction accuracy, and it can be implemented using standard software tools, making it accessible for field engineers and technicians without requiring prior knowledge of computer programming, (5) uncertainty in estimating d_sm is quantitatively evaluated using these literature and novel formulas, and (6) considering 74 formulas from the synthesis of extensive literature gathered and compared in a single study for application according to the need for practical purposes. This enhances the practical applicability and scalability of the developed models in real-world scenarios.

2. Analytical Framework

2.1. Data Description

During the initial stage of data collection, a wide range of more than 40 literature sources are identified. A large collection of more than 1200 unique data samples is obtained during the preliminary stage. Studies with missing data for key parameters were excluded, along with live bed data, which will make significant bed changes [14], focusing only on complete datasets to calculate the dimensionless parameters needed to compute maximum d_s. After such a rigorous screening process focused on important factors influencing scour development, this extremely large dataset is reduced to 768 data samples (#) (Table S4). Here, #: is the number of data points; W: flume width; B: pier diameter; t: time; Δ+1: relative density of sand; H/B: flow shallowness; V/V_c: flow intensity; B/D₅₀: sediment coarseness; Vt/BΔ^0.5: dimensionless time; B/W: constriction ratio. Values are given as minimum/maximum. The key parameter ranges used in the study are assessed overall such as the V/V_c ranged from 0.34 to 1.20, H/B from 0.12 to 21.05, σ from 1.00 to 4.55, B/D₅₀ from 2.25 to 1386.36, Vt/BΔ^0.5 from 7.99×10³ to 3.63×10⁷, B/W from 0.01 to 0.50, F from 0.07 to 1.05, and R from 4.00×10³ to 9.40×10⁵.

2.2. Relationship Between Parameters

In this study, the d_sm for the isolated circular pier is determined using the following set of functional parameters given in Equation (1): -

$$d_{sm}=f\left(H,B,V,V_c,D_{50},{\rho }_f,{\rho }_s,\sigma ,t,g,W,\mu \right)$$

(1)

Here, ρ_f: water density; ρ_s: sediment density; g: gravitational acceleration; µ: dynamic viscosity. The Buckingham π theorem simplifies physical relationships by identifying dimensionless features, reducing variables, and allowing modifications like inverting ratios and converting variables. Adding the relative density parameter (Δ) to non-dimensional time terms is unnecessary, making Equation (1) simplified to Equation (2).

$${\displaystyle \frac{d_{sm}}{B}}=f\left({\displaystyle \frac{H}{B}},{\displaystyle \frac{V}{V_c}},{\displaystyle \frac{B}{D_{50}}},{\displaystyle \frac{{\rho }_{f}VH}{\mu }}=\mathrm{R},{\displaystyle \frac{V}{\sqrt{gH}}}=\mathrm{F},\sigma ,{\displaystyle \frac{Vt}{B{\Delta }^{0.5}}},{\displaystyle \frac{B}{W}}\right)$$

(2)

In this study, the normalized scour depth (d_sm/B) has been used as the predictor, recognizing the inherent flaws and uncertainties associated with equilibrium methods and temporal analysis based on time to attain equilibrium. Thereby, here to counter the uncertainty, the dimensionless time is considered as tV/BΔ^0.5. Existing literatures focused on d_sm under clearwater equilibrium, often neglecting time dependency (Table S5).

2.3. Data Pre-Processing Based on Statistical Measures

During the data pre-processing phase, the original data is transformed into seven independent dimensionless input variables: H/B, V/V_c, B/D₅₀, tV/BΔ^0.5, B/W, F, and R. Issues related to the data pre-processing are given in Section S3 to homogenize the data. The d_sm/B serves as the response variable. Descriptive statistical measures for all dimensionless variables presented in Table S6, alongside the dataset in Table 1, provide a comprehensive overview. The Shields diagram confirms that the selected data are in clearwater condition (Figure S1). The variable tV/BΔ^0.5 exhibits substantial skewness (7.06) and high kurtosis (83.69), indicating deviation from the normal distribution and potential outliers.

The data used in Figure 1 are compiled from previously published experimental studies (Table 1) and supplemented with data from the present laboratory experiment. The correlation matrix illustrates relationships between independent and dependent variables. There is a significant correlation of CC = -0.45 between d_sm/B and B/W, CC = 0.54 between H/B and R, and CC = -0.43 between V/V_c and F. This indicates the importance of H/B, V/V_c, and B/W in predicting d_sm/B. While other predictors display intermediate dependencies, these correlations are based on linear interactions.

2.4. Continuing to Perform Experiments

The research, conducted at Jadavpur University in Kolkata, India, focuses on clearwater scouring. The setup (Figure S3) uses a nonuniform bed material with a 3 m long and 0.81m wide sand bed, aligning mono-piers accurately. Perspex sheets constructed circular piers measuring 5, 7, 9, and 11 cm in diameter were used. Through sieve analysis, sand bed particle sizes D₅₀ = 0.804 mm, σ = 1.75, and (Δ+1) = ρ_s/ρ_f = 2.63 were determined. The experiments lasted 60-78 hours. Flow patterns around the pier and all scour beds post-experiment are analyzed (Figure S4). The methodology flowchart (Figure 2) divides the total data, with 80% used for calibration and 20% for validation.

The bed shear stress (τ_o) =γ_fR_hsinα is derived keeping bed slope (S_s = tanα) constant at 1:2400, R_h: hydraulic radius. Critical bed shear stress (τ_oc) is determined using empirical formulas for the Shields curve, with the Θ_c (= τ_oc/Δρ_fgD₅₀). The critical shear velocity (V_⁎c) is calculated as 0.01966 m/s. Therefore, V_c is calculated from V_c/V_⁎c = 5.75log(H/2D₅₀)+6 where V_⁎c, H, and D₅₀ are already known. The experiments are conducted under clear water scouring conditions (Figure S5). Details of the eight experimental runs and ranges of parameters are set (Table S7) along with the criteria for selecting parameters from the data and literatures (Table S8).

2.5. Machine Learning Models

Gene Expression Programming (GEP) is an evolutionary algorithm that mimics natural selection to evolve computer programs or mathematical functions. It represents potential solutions as linear chromosomes, which are decoded into nonlinear expression trees. The algorithm starts by generating an initial population of random chromosomes. Each chromosome is then expressed as an expression tree, a structure that represents a mathematical function or program. These trees are evaluated based on their performance using a fitness function [94,95]. The fitness of each individual determines their chances of being selected for reproduction. Genetic operations, such as mutation, recombination, transposition, and cloning, are applied to create new offspring. These offspring replace the least fit individuals in the population, forming the next generation. The process repeats until a termination criterion is met, like reaching a maximum number of iterations or achieving a desired fitness level. The best individual from each generation is tracked, with the final best individual representing the solution (Figure 3a). GEP’s fixed-length chromosome structure makes genetic operations easier, while its tree-like expression structure enables the development of complex functions. By combining aspects of genetic algorithms and genetic programming, GEP offers an effective and adaptable method for solving problems like d_s estimation [23,96,97] and other fields [98].

Artificial Neural Network (ANN) model used in this study is a feedforward neural network (FFNN) trained with the Levenberg-Marquardt (LM) backpropagation algorithm. ANNs are widely used for regression and classification tasks due to their ability to model complex, nonlinear relationships between input and output variables [64,99,100,101]. In an FFNN, data moves in one direction, from input to hidden layers, and finally to the output layer, without any cycles or feedback loops [102]. The problem is optimized using a structured hyperparameter tuning approach on an FFNN to achieve accurate predictions. The architecture for the single hidden layer FFNN is given below (Figure 3b).

Key parameters are tuned, including hidden layer size, training epochs, learning rate, activation functions, and training algorithm. By systematically testing different configurations, the model is optimized for low mean squared error (MSE). The LM algorithm is chosen for its efficiency in regression, providing fast convergence with a relatively small dataset. Varying hidden layers, epochs, and learning rates balanced model complexity and stability, allowing the network to learn complex patterns without overfitting. Testing both hyperbolic tangent sigmoid and logarithmic sigmoid activation functions in hidden layers adds flexibility in handling nonlinearity, while pure linear ensures continuous output for regression. Performance is evaluated using MSE across all configurations, tracking the best-performing setup. The final model is then validated on new experimental data, confirming its generalizability.

2.6. Model Evaluation Criteria

The model evaluation criteria, including coefficient of determination (CD), Nash-Sutcliffe efficiency (NSE), mean bias error (MBE), and normalized root mean squared error (RMSE) are employed to assess model performance from different angles: CD checks output alignment with observations, NSE normalizes the statistic to measure relative variances, MBE offers accurate prediction by combining precision and normalization, and RMSE quantifies deviation between observed and measured data. These metrics collectively evaluate efficiency, bias, inaccuracy, and variability, demonstrating their relevance.

3. Results

In this study, an advanced semi-empirical formula for d_sm prediction is developed, providing improved accuracy and adaptability to various hydraulic conditions. A novel GEP model and a robust ANN model are introduced to enhance d_sm prediction accuracy by capturing nonlinear data relationships. The GEP-based formula offers a prediction formula and an interpretable approach, while the ANN model functions as a high-performance black-box solution.

3.1. Semi-Empirical Approach

A semi-empirical formula was developed based on dimensional analysis, where the dimensionless scour depth was modeled as a product of seven key non-dimensional parameters. Each parameter influence was represented by a specific functional form (e.g., sigmoid, Bragg, exponential, logarithmic) chosen based on observed trends. The functions were optimized using nonlinear least squares with the Gauss-Newton method. A detailed explanation of the formulation, parameter fitting process, and optimization steps is provided in Section S4.

The dependence of d_s on seven functional parameters is analyzed. The circle represents the calibration data, and the triangle represents the validation data (Figure 4). The diamond shape points represent the parameter fitted using parameter functions.

Flow intensity (V/V_c): Some researchers found that d_s in uniform sediment increases almost linearly with flow intensity under clearwater conditions, whereas Sheppard and Miller [6] suggested a logarithmic trend. The d_sm reaches the threshold velocity (V/V_c = 1). At this point, the d_s in uniform sediment initially decreases, then rises to a second peak with minor fluctuations [41]. Equation (3) employs two-parameter sigmoid functions in this study to get the best fit of the parameter function of V/V_c, using the Gauss-Newton algorithm for nonlinear regression analysis (Figure 4a). The dependence of d_sm is attained when V/V_c > 0.85.

Flow shallowness (H/B): Various studies [41,103] indicate that when the horseshoe vortex is affected by the creation of the surface roller at the leading edge of the pier, H/B affects d_s as the two vortices rotate in opposite directions. The d_s is therefore independent of H when there is no conflict. Local scour is believed to occur at shallow flow in such instances. As H decreases, the surface roller becomes more prominent, rendering the base vortex less capable of entraining sediment. As a result, with shallower flows, the d_s is decreased. Local scour becomes independent of H in very shallow flows. In such instances, the local scour is considered to occur at a wide pier. The d_sm/B increases with the increase in H/B and yields a maximum at H/B = 4. After this, the parameter function for H/B no longer increases its factor value, which means it has no effect on computing d_sm given in Figure 4(b). However, some researchers found maximum d_sm/B at varying H/B ratios, like H/B =3 [104], H/B = 10 given by Sheppard and Miller [6], etc. This increasing trend has already been reported in terms of power law [45,105]. In the present study, two-parameter sigmoid functions are applied to get the best fit of the parameter function H/B given in Equation (4), using the Gauss-Newton algorithm for nonlinear regression analysis (Figure 4b). The analysis reveals that d_sm is independent when H/B > 4.

Sediment gradation (σ): Non-uniform sediment can result in a lower d_s and an exponential declining tendency with increasing σ [19]. The highest value produced for homogenous sediment is about σ=1-1.3, with a declining tendency after σ =1.5. The effect of σ on scour was also used in predicting d_s mentioned earlier by several researchers. The effect of σ on the dependence of d_sm diminishes when σ>2.5, and the dependence of d_sm on σ is shown using Bragg’s function given in Equation (5) (curve with maxima and minima) (Figure 4c).

Sediment coarseness (B/D₅₀): Variation in d_s remains unaffected by B/D₅₀ in the presence of large-scale turbulence for homogeneous sediments, unless D₅₀ is relatively large. Ettema [68] observed that when B/D₅₀ is low (<50), individual grains are large relative to the groove excavated by the downflow, and erosion is hampered because the sediment bed dissipates some energy of the downflow. When B/D₅₀<8, individual grains are so large that scour is primarily due to erosion along the pier sides. Sheppard et al. [79] showed considerable d_s reductions for increasing B/D₅₀ at much larger values of B/D₅₀, which differ from the value found by other researchers. Pandey et al. [21] preferred to use a single straight line to illustrate the effect of d_sm/B with B/D₅₀ continuously increasing with B/D₅₀ upto B/D₅₀=650, whereas Franzetti et al. [14] used an exponential function to show an increasing trend d_sm/B with B/D₅₀ upto B/D₅₀=60 followed by a decreasing trend given in Equation (6). This study demonstrates an exponential relationship to illustrate the relation between d_sm/B and B/D₅₀ (Figure 4d). The figure shows an upward trend in the parameter function when B/D₅₀ <55, transitioning to a downward trend when B/D₅₀ ≥ 55. Beyond B/D₅₀ >760, the effect on the relation between d_sm/B and B/D₅₀ diminishes, suggesting a limiting condition where sediment size no longer governs the scour process significantly.

Time scale (tV/BΔ^0.5): Some researchers employed the time scale as tV/B to compute d_st as a logarithmic function [106]. The particular impact of Δ is not evaluated in the current analysis. Instead, it is included in the time scale (=tV/BΔ^0.5), serving as a controlling factor for temporal scour. Some researchers employed exponential trends in the parameter function of tV/BΔ^0.5 [13]. Therefore, a relationship in the dependence between d_sm/B and tV/BΔ^0.5 is given using an exponential correlation (Figure 4e) and Equation (7).

Constriction ratio (B/W): Blockages from side walls in flume tests have minimal effects on scour when the channel is sufficiently large, with a B/W < 10% [70]. Conversely, for B/W >10%, the side wall has a significant influence on flow characteristics and scour. On B/W impact, Sheppard et al. [79] researched for B/W of 15%, while Aksoy et al. [89] focused on B/W > 10%. An exponential relationship of dependence between d_sm/B and B/W is shown (Figure 4f) and in Equation (8). In this study, it is observed that the impact of the B/W on the dependence of d_sm/B is minimal.

Flow Froude number (F): Different researchers applied F to estimate d_s considering the power law [107]. While others employed either densimetric Froude number (F_d) as a scour-affecting parameter [34], or critical Froude number (F_c) as the governing parameter [108]. In this study, the two-parameter sigmoid function is used using the Gauss-Newton algorithm for nonlinear regression analysis to obtain the best fit of the parameter function of F given in Equation 9 (Figure 4g). The figure confirms the dependence of d_st/B on F attains when F = 0.3, and that effect tends to diminish when F> 0.3.

The impact of Δ on time is considered following Franzetti et al. [14]. While the influence of B/W is minor, the role of W, mentioned in Sheppard et al. [79], remains important. After evaluating, it is found that R has little impact on d_sm calculation. Hence, R is omitted from d_sm predictor function list.

The impact of F is considerable within lower value ranges. As it increases, its influence diminishes significantly in calculating d_sm. The scour function is defined by η = 2.963, with its maximum value occurring when all parameter function values reach 1. Essentially, η indicates d_sm at points of maximum parameter function values (Equations (3)–(9)).

$$f_1\left({\displaystyle \frac{V}{V_c}}\right)=1-\exp \left[-\exp \left\{-3.35949+6.03486\left({\displaystyle \frac{V}{V_c}}\right)\right\}\right]$$

(3)

$$f_2\left({\displaystyle \frac{H}{B}}\right)=1-\exp \left[-\exp \left\{-0.27283+0.55173\left({\displaystyle \frac{H}{B}}\right)\right\}\right]$$

(4)

$$f_3\left(\sigma \right)=0.35379+\left(0.99199-0.35379\right)\left[\exp \left\{-2.12547{\left(\sigma -1.13131\right)}^2\right\}\right]$$

(5)

$$f_4\left({\displaystyle \frac{B}{D_{50}}}\right)=0.8491\left[{\left({\displaystyle \frac{B}{D_{50}}}\right)}^{1.8158}\exp \left\{-2.9905{\left({\displaystyle \frac{B}{D_{50}}}\right)}^{0.235}\right\}+0.5112\right]$$

(6)

$$f_5\left({\displaystyle \frac{tV}{B{\Delta }^{0.5}}}\right)=1-\exp \left\{-0.09000{\left({\displaystyle \frac{tV}{B{\Delta }^{0.5}}}\right)}^{0.20293}\right\}$$

(7)

$$f_6\left({\displaystyle \frac{B}{W}}\right)=\exp \left\{-1.03251{\left({\displaystyle \frac{B}{W}}\right)}^{2.50001}\right\}$$

(8)

$$f_7\left(\mathrm{F}\right)=1-\exp \left[-\exp \left\{-0.362465+8.20785\left(\mathrm{F}\right)\right\}\right]$$

(9)

The new scour predictor formula is represented in Equation (10) as the product of the effect of the individual parameters [109] and optimizing each parameter function by combining Equations (3-9).

$${\displaystyle \frac{d_{sm}}{B}}=\eta f_1\left({\displaystyle \frac{V}{V_c}}\right)f_2\left({\displaystyle \frac{H}{B}}\right)f_3\left(\sigma \right)f_4\left({\displaystyle \frac{B}{D_{50}}}\right)f_5\left({\displaystyle \frac{Vt}{B{\Delta }^{0.5}}}\right)f_6\left({\displaystyle \frac{B}{W}}\right)f_7\left(\mathrm{F}\right)$$

(10)

3.2. Machine Learning Approach

3.2.1. Gene Expression Programming

The prediction of d_s around bridge piers was developed using the GEP model applied to the collected big dataset and the present experimental dataset. To check the accuracy of GEP models in the deployment phase, in the MS Excel workbook, which is divided into a training subset (614 observations, 79.95% of the total) and a testing subset (154 observations, 20.05% of the total), to ensure model validation. The RMSE is selected as the fitness function for evaluating the model’s accuracy, a choice supported by prior studies [23]. A population size of 30–100 is used for optimal performance, as suggested by Ferreira [94], and RMSE is calculated for each chromosome to determine its fitness. The formulation process involves five primary steps. First, the initial population is established, followed by fitness evaluation based on RMSE. Then, a selection of functions and terminals is made, drawing on both mathematical operators (Table 2) and relevant dimensional variables to capture influential parameters in d_s prediction. The architecture of each chromosome is carefully defined. In the final steps, genetic operators such as gene transposition, inversion, and recombination are applied to optimize the solution.

The GEP model ultimately generates an explicit equation for d_sm (Equations (20-21)) as a function of dimensionless parameters given in Equation (2). This equation is represented through three expression trees (ETs), each linked by an addition function to account for cumulative contributions to d_s (Figure 5). The proposed d_s prediction equation combines dimensional factors crucial to scour processes, providing a robust predictive tool supported by both training and test data. The final GEP-based equation serves as a reliable predictor of d_s, with its parameters and architecture informed by extensive experimental data and validated through rigorous testing of literature and with additional experimental data.

$${\displaystyle \frac{d_{sm}}{B}}=\left\{\begin{array}{l}\max \left[d_6,\log \left\{\tanh \left(d_1\right)+{\displaystyle \frac{d_0}{d_2}}\right\}+\min \left(d_0+G1C1,{\displaystyle \frac{1}{d_2}}\right)\right]\\ \times {\left[G2C9\ast \tanh \left\{{\displaystyle \frac{\frac{d_3}{G2C5}-\mathrm{atan}\left(d_1\right)}{2}}+\min \left({\displaystyle \frac{d_2}{d_1}},d_1\right)\right\}\right]}^2\\ \times {\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{1}{\max \left(d_0,G3C2\right)}}-\log \left({\displaystyle \frac{d_2+d_0}{2}}\right)+{\displaystyle \frac{\log \left(d_6\ast d_4\right)+\left(G3C9\ast d_5\right)}{2}}\right]\end{array}\right.$$

(11)

$${\displaystyle \frac{d_{sm}}{B}}=\left\{\begin{array}{l}\max \left[\mathrm{F},\log \left\{\tanh \left({\displaystyle \frac{H}{B}}\right)+{\displaystyle \frac{\left(\frac{V}{V_c}\right)}{\sigma }}\right\}+\min \left({\displaystyle \frac{V}{V_c}}+0.03795,{\displaystyle \frac{1}{\sigma }}\right)\right]\\ \times {\left[0.6467\ast \tanh \left\{{\displaystyle \frac{\frac{\left(\frac{B}{D_{50}}\right)}{-5.9929}-\mathrm{atan}\left(\frac{H}{B}\right)}{2}}+\min \left({\displaystyle \frac{\sigma }{\left(\frac{H}{B}\right)}},\left({\displaystyle \frac{H}{B}}\right)\right)\right\}\right]}^2\\ \times {\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{1}{\max \left(\frac{V}{V_c},0.8306\right)}}-\log \left({\displaystyle \frac{\sigma +\frac{V}{V_c}}{2}}\right)+{\displaystyle \frac{\log \left(\mathrm{F}\ast \frac{tV}{B{\Delta }^{0.5}}\right)+\left\{\left(-18.7964\right)\ast \left(\frac{B}{W}\right)\right\}}{2}}\right]\end{array}\right.$$

(2)

3.2.2. Feedforward Neural Network (FFNN)

An FFNN is an ANN where data flows in one direction from the input through hidden layers to the output layer. FFNNs are valued for their simplicity and efficiency in generating direct predictions, as information moves only forward through the network.

Training the FFNN begins with random weights, which are iteratively adjusted to minimize prediction errors through a process called back-propagation. This process fine-tunes weights and biases to improve model accuracy. To train the FFNN in this study, the data is split into three sets: training (70%), validation (15%), and testing (15%). The training set adjusts the model’s weights and biases, the validation set helps avoid overfitting, and the testing set evaluates the model’s performance. Three optimization algorithms, Levenberg-Marquardt (LM), are used to tune the FFNN model. The LM algorithm, available in MATLAB, is implemented without further customization, while the grid search structural approach is specifically explained and applied for parameter optimization. The learning of FFNN and training progress are shown (Figure 6).

Figure (6a) illustrates the training optimization process of FFNN in terms of the accuracy measure of MSE and Epochs, 28 is found to be the optimal value. A decreasing gradient signifies that the model is effectively minimizing the error between its predictions and the actual values. The middle plot shows the learning rate, often denoted as “mu” (Figure 6b). This parameter controls the step size taken during each update to the parameters. A relatively stable learning rate suggests a well-tuned optimization process. The bottom plot visualizes the number of validation checks performed at each epoch. These checks help assess the model’s performance on unseen data and prevent overfitting. The increasing trend in validation checks towards the end of training is a common practice to ensure the model’s generalizability. Overall, the graph indicates a successful training process where the model gradually improved, likely achieving satisfactory performance on both the training and validation datasets.

Figure 7 shows four scatter plots, each representing the relationship between the predicted and observed values of a model. The top left plot shows the training data, where the model’s predictions (blue line) closely follow the actual data points with CD=0.8338. The top right plot shows the validation data, where the model’s predictions (green line) are slightly less accurate but still reasonably close to the actual data points with good CD=0.7727.

The bottom left plot shows the test data CD=0.7901, where the model’s predictions (red line) are less accurate than the training and validation data. The bottom right plot shows all the data combined, where the model’s predictions (cyan line) are a compromise between the training, validation, and test data (CD=0.8182). Overall, the figure suggests that the model is performing reasonably well.

Figure 8 shows a histogram of errors for FFNN. The histogram is divided into 20 bins, each representing a range of error values. The bars in the histogram show the number of instances (data points) that fall within each error range. The blue bars represent the training data, the green bars represent the validation data, and the red bars represent the test data. The orange line indicates zero error, which means the model’s prediction was exactly correct. It can be seen that for the training data, most errors are clustered around zero, indicating good performance. However, the validation and test data show a wider spread of errors, suggesting that the model might be overfitting to the training data.

3.3. Evaluation of Models/Formulas

Evaluation matrices including CD, NSE, MBE, RMSE, and percentage of data within the desirable error limit (P_in) for training/calibration, validation, and experimental datasets are shown in Table 2 for the semi-empirical model, GEP, and FFNN. The highest-performing models in each dataset category highlight the strengths of different modeling approaches. In the training/calibration, validation, and testing phase, FFNN demonstrates superior predictive accuracy. For the present experimental data, the GEP model leads, showing strong performance. This result suggests that FFNN excels in controlling all scenarios except for the additional experimental data. But in the case of P_in, the empirical and FNN models for experimental data perform best. Other than all cases, FFNN performs best in terms of P_in. It shows the ranking of formulas based on the percentage of data falling within the error zone of ±20% of the best agreement line using the normalized value. Figure 9 indicates that the new formula is better than all 74 literature formulas with respect to most of the datasets from 40 literatures (Table S9).

Based on the ranking using P_in, and comparison d_sm/B between measured and experimental d_sm values from the top three literatures along with the present formulas (Figure 9). It is seen that using FFNN, most of the data are closer to the perfect agreement line (PAL), indicating a better model among others, including the literature formula (Figure 9a). However, GEP performs and generalizes better than other formulas for the present experimental datasets. The semi-empirical model of the present study performs better than other literature formulas (Figure 9c), such as the formula by Nandi and Das [15], Franzetti et al. [14], and Lança et al. [13]. The overestimation of d_sm/B can be seen for Lança et al. [13] (Figure 9f).

A comparison of different models or methods using a combined box and violin plot, called a box-violin plot, is shown below (Figure 10). The observed data has a certain median and spread, which other models ideally aim to match closely, as this would suggest high accuracy in their predictions. The median of FFNN aligns well with the observed median, and the spread is similar, suggesting that FFNN is effective. The median GEP is close to the observed data, and the spread bit thicker than FFNN and observed values suggest it is less accurate than FFNN. While the median and spread of Lança et al. [13] are not similar to the observed significant differences in the distribution or median position, indicating limitations in the model’s applicability to newer or different datasets.

The FFNN and GEP models give results that are very close to the observed values, with less spread and better consistency. The Semi-Empirical model also shows good performance, better than the older models. In comparison, the models by Franzetti et al. [14] and Lança et al. [13] show more variation compared to observed data. Overall, FFNN performs the best, followed by GEP and the Semi-Empirical model, showing that the models used in this study are more reliable than the older ones (Figure 10).

The Taylor diagram also shows that the FFNN models give better results as it closer to the observed values and concentrate on the higher correlation line, while Lança et al. [13] perform worse as the point is situated far away from the observed value and concentrated towards the lower correlation line. However, the performance of GEP is slightly better than the semi-empirical formula shown (Figure 11).

3.4. Overall Discussion

The predicted values of d_sm/B for additional experimental data using FFNN, GEP, and semi-empirical models align well with observed values. Nandi and Das’s [15] formula also shows strong agreement between calculated and measured results, ranking in 4th place in terms of overall data points inside the 20% error band. The majority of literature datasets show a strong agreement, highlighting consistency between the present predicted d_sm/B with observed d_sm/B across different literature. The comparison of formulae from all literature sources, along with plots, is detailed (Figure S6).

It is observed that the present formula (Semi-Empirical) has 2.5% more influence based on CD, 4.6% better accuracy based on NSE, 28.3% better according to MBE, and 4.9% more accurate with respect to RMSE in predicting d_sm compared to the best literature formula (Table S9).

It is observed that, comparing the overall data, GEP has 5.1% more influence based on CD, 4.4% better accuracy based on NSE, 83.3% better according to MBE, and 9% more accurate with respect to RMSE in predicting d_sm compared to the best empirical formula (present study).

It is also observed that the present formula (FFNN) has 10.1 % more influence based on CD, 10.9% better accuracy based on NSE, 28.3% better according to MBE, and 6.2% more accurate with respect to RMSE in predicting d_sm compared to the GEP. The statistical accuracy of the remaining formulas is found to be significantly less than the newly proposed formula (Table S9).

In the U.S., the commonly used HEC-18 formula was initially proposed by Richardson and Davis [1], and then revised by Arneson et al. [7] based on experiments. A comparison in Table S9 shows that, for extensive datasets, the HEC-18 formula outperforms the revised version. Gao et al. [36] developed a formula using field data, which has been used in Chinese highway and railway engineering for over two decades. However, the statistical influence parameters of this formula have been shown to be inadequate, with a P_in of only 35.18%.

The Indian practice of estimating design d_s, known as Lacey’s method, was developed by Lacey [32] and Inglis [33] based on observations made on canals in India and Pakistan. This method is commonly used to estimate d_s around bridge piers in alluvial rivers and is approved for design by Indian Railways and IRC-78 [29]. However, the statistical parameters for IRC-78 [29] do not show better results with a low P_in of 26.16%. Vijayasree and Eldho [35] presented a modified IRC formula, which provides better predictions for different shapes. They suggested the modified IRC formula for bridge foundation design under Indian conditions. The statistical parameters for that formula do not work well with a very low P_in of just 6.701% (Table S9). Therefore, in comparing IRC-78 [29] and modified IRC using the extensive data in this study, IRC-78 [29] significantly outperforms modified IRC.

However, while each formula has advantages, it also has its drawbacks. HEC-18 may lack accuracy in some scenarios; the Gao et al. [36] formula may not fully capture all the variables that affect d_s. Although the Lacey method is widely used in India, it cannot accurately estimate the characteristics of alluvial rivers. Furthermore, the modified IRC formula, although promising, may require further validation to confirm its applicability in diverse conditions.

Increasing the quantity of experimental data may increase the power of the analysis and provide a greater understanding of the parameter being considered. The proposed model has not yet been applied to real-world case predictions, as it is developed purely based on controlled, clear water laboratory experimental data, which is the main limitation of this study. To make the model suitable for practical field applications in the future, it would need to be retrained the models using observed real-world data. In addition, proper consideration of scaling laws and field conditions would be necessary to ensure the accuracy of models and reliability outside the laboratory environment. Field data may be taken into consideration to assess the applicability of current results in real-world situations, in order to offer insightful information on the scour parameters. Other advanced machine learning models with different search optimization techniques may give better results.

3.5. Uncertainty Analysis

This section compares the newly proposed formula (GEP and semi-empirical)/models (FFNN) with literature formulas to quantitatively evaluate uncertainty in estimating d_sm. This study uses an uncertainty analysis database containing 768 experimental measurements taken from literatures to calculate individual prediction errors. The calculated mean and SD of these errors indicate that the observed values are overestimated for positive means and underestimated when negative means appear. At 5% significance level, this study (FFNN) exhibits the lowest prediction uncertainty band (0.037) among all formulas considered (Table S10). These results offer valuable insights into the reliability and accuracy of the new prediction formula, which improves the understanding of d_sm prediction as evidenced by the uncertainty analysis.

4. Conclusions

Scouring around bridge piers can weaken the structure. If not resolved, scour-related problems may cause bridge accidents, impacting public safety, transportation, and the economy. In this study, novel formulas (GEP and semi-empirical)/models (FFNN) have been developed to assess maximal scour depth (d_sm), incorporating methodologies applied to datasets from 40 different literature sources. In computing the normalized maximal scour depth (d_sm/B), the best input parameters were found as flow intensity (V/V_c,), flow shallowness (H/B), sediment gradation (σ), sediment coarseness (B/D₅₀), dimensionless time (tV/BΔ^0.5), constriction ratio (B/W) and Froude number (F). The performance of new formulas/models was evaluated using statistical indices: CC, NSE, MBE, and RMSE. The accuracy of novel formulas was tested against more than 74 empirical literature formulas.

• o Two additional parameters, B/W and F, significantly influence the d_sm prediction within ranges from 0.1 to 0.5 and 0.08 to 0.28, respectively.
• o The FFNN achieved the highest performance in both testing (CD = 0.790, NSE = 0.783, MBE = -0.039, and RMSE = 0.289) and validation (CD = 0.773, NSE = 0.767, MBE = -0.040, and RMSE = 0.266. In comparison, the semi-empirical formula and the GEP model.
• o Newly developed FFNN models provide 18.6% more accuracy in terms of CD and 15% in terms of P_in (%) predictions compared to the most accurate literature formula.
• o The present study (GEP) formula provides 7.7% better CD compared to the existing best empirical formula. It gave steady and reliable results, with very little difference in performance measures across calibration, validation. This consistency shows its reliability across various datasets.
• o The present study (Semi-empirical) formula outperformed all the existing literature formulas.
• o Furthermore, the uncertainty analysis performed shows that FFNN gives the best results with the least margin of error (0.018) and the narrowest computation uncertainty band (0.037).

The limitation is that the predictive abilities of Semi-empirical, GEP, and ANN models may not work well for new field data or live bed conditions that are very different from the present calibration data, which relies solely on laboratory-based data. The different advanced machine learning models may be used for better predictive accuracy. The field data may be added to capture the complexity of real-world or natural flow conditions and generalise the model for practical implementation in the future. The advanced optimization technique may improve the predictive performance in the future. More research is needed to expand the predictive range of these models, making them reliable for both lab and field applications, beyond just clear water data.