Estimation and interpretation of equilibrium scour depth around circular bridge piers by using optimized XGBoost and SHAP

Abstract

Most bridge failures result from scouring around bridge piers, resulting in economic losses and risks to public safety. The conventional equations for predicting the depth of scour at bridge piers have several limitations: (1) They mainly use regression-based techniques that cannot robustly capture the nonlinear relationship between the scour depth and its effective variables; (2) they are applicable only to a narrow range of variability of data; and (3) they are typically calibrated using laboratory data rather than field measurements and thus cannot simulate the prototype environment. To overcome these limitations, in this study, three novel hybrid machine learning methods: particle swarm optimization - extreme gradient boosting (PSO – XGBoost), red fox optimization - XGBoost (RFO – XGBoost), and relativistic particle swarm optimization - XGBoost (RPSO – XGBoost) are applied to estimate the scour depth around circular bridge piers, and their effectiveness is validated using three statistical metrics, i.e. the mean absolute error (MAE), root mean square error (RMSE), and coefficient of determination (R²). RPSO – XGBoost generates the best results for both dimensional and dimensionless data. Moreover, the proposed approaches outperform the state-of-the-art techniques. The SHapley Additive exPlanations (SHAP) method is used to assess the relative significance of the contributing factors for predicting the scour depth.

KEYWORDS:

• Scour depth
• XGBoost
• Hyperparameter optimization
• Explainable AI
• Bridge pier
• SHAP

1. Introduction

Scouring around bridge piers is a notable driving factor of bridge collapses (Wardhana & Hadipriono, 2003; Wang et al., 2017). According to the US Department of Transportation, scouring led to the damage of 17 bridges in New York and New England during the flood in 1987 (Richardson & Davis, 2001). The Federal Highway Administration defines scouring as the process in which the streambed sediments are eroded and removed from bridges by the water flowing over the bed (Mueller & Wagner, 2005). Different bed materials undergo various levels of scouring, leading to different scour rates. For example, non-cohesive materials such as sand and gravel may erode over a few hours, but more time is required for cohesive bed materials to be eroded (Richardson & Davis, 2001).

The scouring process around bridge piers is composed of two mechanisms: three-dimensional flow separation and sediment transport, which render it complex (Sreedhara et al., 2019). The water rushing around the pier experiences three-dimensional separation, which generates pressure resulting in a downward flow upstream of the pier. This downward flow creates a vortex system upstream of the pier, termed a horseshoe vortex owing to its appearance being similar to that of a horseshoe from the top view (Pandey et al., 2020a). The shear stress increases due to the vortex system and corresponding downflow around piers, leading to accelerated sediment transport, which finally creates a scour hole near the pier. After the scour hole is fully developed and reaches its maximum size, the flow pattern is altered, and the shear stress and sediment transport diminish (Pandey et al., 2020b).

The safe design of bridge foundations is critical for preventing bridge failures. The scour depth must be precisely estimated to ensure the economic and secure design of bridges. Hence, many regression-based formulas have been established to determine the depth of equilibrium scour around bridge piers (Laursen & Toch, 1956; Shen et al., 1969; Hancu, 1971; Jain & Fischer, 1979; Kothyari et al., 1992; Lee & Sturm, 2009; Melville & Coleman, 2000; Khan et al., 2017; Pandey et al., 2018). Notably, the existing equations suffer from the following limitations: (1) They are mainly based on limited experimental data and are thus applicable to only a narrow range of variability of inputs (e.g. flow depth, flow velocity, pier diameter, and sediment size) and output (scour depth) for which they are calibrated; (2) they cannot robustly capture the highly nonlinear relationship between scour depth and its determining factors; and (3) they are typically calibrated based on laboratory data and thus cannot accurately capture the prototype environment, resulting in overestimation of the scour depth (Bateni et al., 2007a).

In a departure from regression-based approaches, machine learning (ML) techniques have been used to predict scour depth. Sharafi et al. (2016) showed that the support vector machine (SVM) outperforms artificial neural networks (ANN), adaptive neuro-fuzzy inference systems (ANFIS), and nonlinear regression methods in estimating the scour depth. Ebtehaj et al. (2017) compared the scour depth prediction performance of the self-adaptive extreme learning machine with those of regression-based, ANN, and SVM techniques. Choi et al. (2017) used ANFIS to determine the scour depth around circular bridge piers. While their approach outperformed the conventional equations, the small number of data points in their study restricts the generalizability of their proposed method. Genetic programming and ANFIS were proposed as alternatives to conventional methods owing to their promising performance in estimating the scour depth near piers (Abd El-Hady Rady, 2020). Shamshirband et al. (2020) developed several PSO-based equations to predict the scour depth around circular bridge piers using different input variables.

Dang et al. (2021) optimized ANN by PSO and firefly algorithms to enhance its efficiency in scour depth estimation. However, the performances of their proposed methods were not compared with other studies, limiting the ability to assess their effectiveness. Sreedhara et al. (2021) utilized gradient tree boosting (GTB) to estimate the temporal scour depth around circular, rectangular, sharp-nosed, and round-nosed bridge piers. Their proposed ML model was calibrated based on a limited number of laboratory data points and was not evaluated against previous studies. Chou and Nguyen (2022) introduced the Metaheuristic-Optimized Stacking System (MOSS) to estimate the scour depth around bridge piers. By combining a metaheuristic algorithm with multiple ML methods, they showed the superiority of MOSS over single ML models, conventional regression-based methods, and mathematical approaches. Despite the strength of their study, it used a small dataset, which hinders its generalization. Also, none of the abovementioned studies used explainable artificial intelligence (AI) to show the relative contribution of input features to scour depth estimations.

This study overcomes the drawbacks of the existing studies by (1) using 841 experimental data points from 35 field and laboratory studies in clear water conditions, thereby covering a wide range of variability of the scour depth and other relevant variables; (2) applying three novel hybrid ML techniques: particle swarm optimization–extreme gradient boosting (PSO–XGBoost), red fox optimization–XGBoost (RFO–XGBoost), and relativistic particle swarm optimization–XGBoost (RPSO–XGBoost), given that PSO, RFO, and RPSO can tune the hyperparameters of XGBoost to improve its performance; and (3) interpreting the proposed ML models using the SHapley Additive exPlanations (SHAP) method. In particular, SHAP represents an intuitive and effective tool for exploring the contribution of each input feature to the output of an ML model.

XGBoost is a robust and scalable ML method that was proposed by Chen et al. (2015). This framework operates on the basis of a gradient boosting (GB) tree for gradient enrichment. XGBoost has demonstrated robust performances in regression problems, in terms of speed, memory usage, scalability, and hardware (Lucca et al., 2021; Qiu et al., 2021), in applications such as detection in water delivery systems (Wu et al., 2022), evaluation of flash flood risk (Ma et al., 2021), groundwater level estimation (Osman et al., 2021), and flood forecasting (Venkatesan & Mahindrakar, 2019). In general, XGBoost can model a variety of water resource problems, and its performance can be improved by tuning its hyperparameters (Ni et al., 2020; Yu et al., 2020; Shi et al., 2021; Nguyen et al., 2021; Demir & Sahin, 2023). Although this tuning can be performed through trial and error, such approaches are extremely time-consuming and cumbersome.

PSO is a popular swarm-based metaheuristic optimization technique that draws inspiration from animal foraging behaviours (Kennedy & Eberhart, 1995). RFO, as one of the newest metaheuristic optimization tools, was introduced by Połap and Woźniak (2021). Its concept is based on the predation tricks of red foxes (Cui et al., 2022; Natarajan et al., 2022). RPSO, introduced by Roder et al. (2020), is a variant of the PSO algorithm that functions based on the concept of relativity. These optimization strategies have been successfully used to adjust hyperparameters of ANN (Alizamir & Sobhanardakani, 2018), XGBoost (Yu et al., 2020; X. Zhang et al., 2020), and random forest platforms (Pham et al., 2020).

In this study, we exploit the PSO, RFO, and RPSO techniques to identify the optimal hyperparameters of XGBoost to improve its performance (Gu et al., 2021; Lucca et al., 2021). To the best of our knowledge, the proposed hybrid approaches have not been previously used to predict the scour depth. Moreover, the performances of the developed hybrid models are compared with those of state-of-the-art methods.

Furthermore, the interpretation of ML methods is crucial to explain how these algorithms obtain optimal predictions. To this end, explainable artificial intelligence (XAI) methods have attracted increasing attention in recent years. Representative XAI techniques include SHAP, accumulated local effects, and partial dependence plots (Ishfaque et al., 2022).

In this study, the relative significance of the different input parameters for scour depth prediction is determined using SHAP values. SHAP, proposed by Lundberg and Lee (2017), quantifies the significance of each feature and explains how ML models make predictions. The predictions of a model can be expressed as the summation of a fixed base value and each feature’s corresponding SHAP values. Unlike the conventional feature importance methods such as feature importance and mean decrease impurity (Mangalathu et al., 2020; Demir & Sahin, 2023), SHAP can determine whether the contribution of each feature is positive or negative. SHAP has been used to interpret several ML methods for dam seepage problems (Ishfaque et al., 2022) and drought forecasting (Dikshit & Pradhan, 2021) owing to its notable advantages such as local and global model interpretation, scalability, robustness, and feasibility for a vast range of problems such as regression, classification, and ranking (Lundberg et al., 2020). To the best of the authors’ knowledge, this study represents the first attempt at performing SHAP analysis to interpret ML methods for scour depth prediction.

The rest of the paper is organized as follows: Section 2 describes the variables affecting the scour depth and data sources. Section 3 discusses the XGBoost, PSO, RFO, RPSO, and SHAP techniques. Section 4 presents the findings, and Section 5 presents the concluding remarks.

2. Data

The scour depth, influenced by flow characteristics, sediment properties, and pier dimensions, is inherently complex due to the intricate interplay of hydraulic forces, sediment dynamics, and morphological changes (Choi et al., 2017). This complexity accentuates the need for extensive understanding and modeling techniques to analyze and predict the scour depth accurately. However, conducting three-dimensional physical-based simulations around bridge piers to accurately represent the scour mechanism is challenging (Dang et al., 2021). Given the physical complexity and the limitations of physical models, the availability of a comprehensive scour dataset becomes crucial. When utilized with ML approaches, such a dataset provides valuable insights into the intricate scour phenomena and enables the rapid generation of accurate scour depth estimates (Dang et al., 2021). Importantly, our study stands out as the first one to employ such a comprehensive dataset to model scour depth around circular bridge piers, further highlighting its significance.

The maximum scour depth, also known as the equilibrium, scour depth, is reached around bridge piers after a period of gradual development. The equilibrium scour depth near circular bridge piers is determined by the fluid flow, bed materials, and pier diameter (Melville & Coleman, 2000): (1) $d_{se}=f\left(\rho ,\mu ,V,Y,g,d_{50},V_c,D\right)$(1) where d_se is the equilibrium scour depth, ρ is the density of water, μ is the dynamic viscosity of water, V is the flow mean velocity, Y is the flow depth, g is the gravitational acceleration, d₅₀ is the median sediment size, V_c is the mean critical velocity, and D is the pier diameter.

Using the Π Buckingham theory and considering ρ, D, and V as repeating parameters, the dimensionless form of Eq. (1) can be expressed as follows: (2) $\frac{d_{se}}{Y}=f\left(\frac{V}{V_c},\frac{V}{\sqrt{gY}},\frac{D}{Y},\frac{d_{50}}{Y},\frac{\rho VD}{\mu }\right)$(2) Here, $\frac{V}{\sqrt{gY}}$ is the Froude number (Fr), and $\frac{\rho VD}{\mu }$ is the Reynolds number (Re).

Reynolds number determines the flow regime around bridge piers, where a higher Reynolds number signifies a more turbulent condition that may enhance the erosive potential and increases scouring (Hu et al., 2022). Therefore, including Re in simulating scouring around bridge piers can be beneficial to account for its underlying effect and ensure a thorough understanding of the scouring process (Tavouktsoglou et al., 2017). However, it is worth mentioning that our results in Section 4.2. show that Re has a marginal impact on the scour depth estimates, and its effect on d_se/ Y is less than the other variables in equation (2).

In this study, 841 field and laboratory data points in clear water condition for circular bridge piers are collected from 35 published and unpublished reports, which are listed in Table 1.

Table 1. Sources for scour depth data used in this study.

Study	Type of data	Number of data points
Chabert and Engeldinger (1956)	Laboratory	87
Shen et al. (1969)	Laboratory	2
Hancu (1971)	Laboratory	1
Ettema (1976)	Laboratory	19
Jain and Fischer (1979)	Laboratory	6
Ettema (1980)	Laboratory	97
Jain (1981)	Laboratory	2
Chiew (1984)	Laboratory	11
Yanmaz and Altinbilek (1991)	Laboratory	33
Kothyari et al. (1992)	Laboratory	92
Dey et al. (1995)	Laboratory	18
Graf (1995)	Laboratory	3
Melville (1997)	Laboratory	17
Melville and Chiew (1999)	Laboratory	27
Oliveto and Hager (2002)	Laboratory	44
Mia and Nago (2003)	Laboratory	5
Sheppard (2004)	Laboratory	12
Raikar and Dey (2005)	Laboratory	20
Ettema et al. (2006)	Laboratory	6
Kumar (2007)	Laboratory	6
Unpublished (Bateni et al., 2007a)	Laboratory	33
Unpublished (Bateni et al., 2007a)	Field	27
Coleman (unpublished) (Sheppard et al., 2010)	Laboratory	6
Jones (unpublished) (Sheppard et al., 2010)	Laboratory	16
Sheppard et al. (2010)	Laboratory	5
Kumar (2011)	Laboratory	3
Lança et al. (2013)	Laboratory	38
Lodhi et al. (2014)	Laboratory	3
López et al. (2014)	Laboratory	29
Aksoy and Eski (2016)	Laboratory	28
Fael et al. (2016)	Laboratory	4
Aksoy et al. (2017)	Laboratory	16
Williams et al. (2017)	Laboratory	22
Pandey (unpublished)	Laboratory	18
Pandey et al. (2020a)	Laboratory	85

Table 2 summarizes the mean, minimum value, maximum value, standard deviation (SD), and coefficient of variation (CV) of the data.

Table 2. Statistical indices of variables.

Variable	Minimum	Maximum	Mean	Standard deviation (SD)	Coefficient of variation (CV)
d50 (mm)	0.220	14.250	2.333	2.619	1.122
D (m)	0.010	1	0.129	0.153	1.194
V (m/s)	0.149	1.270	0.417	0.223	0.535
Vc (m/s)	0.223	1.274	0.525	0.243	0.464
Y (m)	0.020	1.900	0.231	0.215	0.927
dse (m)	0	1.410	0.129	0.124	0.964

3. Models and methods

3.1. XGBoost

XGBoost is a novel variant of the GB technique (Friedman, 2001, 2002). Unlike the GB technique, XGBoost can efficiently manage a considerable amount of data by generating boosted trees and implementing them in parallel (Le et al., 2019; Zhang et al., 2020). Consequently, reliable and fast simulations can be performed for engineering problems with a large number of data points. The complexity of trees is addressed by the variation of a loss function. In other words, the residuals are used to calibrate the former predictions in each iteration through the optimization of a loss function. Furthermore, XGBoost uses a regularization function in the objective function (Zhou et al., 2021).

The process flow of XGBoost can be mathematically illustrated as follows: $\begin{array}{rl}& \text{Consider a dataset with }m\text{ data and }n\text{ feature}:\\ & S=\left(x_i,y_i\right)\left(|S|=m,x_i\in R^n,y_i\in R\right).\end{array}$x_i and y_i are the input and output, respectively. The output of a tree ensemble model can be obtained by adding P functions: (3) ${\hat{y}}_i=\psi \left(x_i\right)=\sum _{p=1}^P{f}_p\left(x_i\right),f_p\in F$(3) where F is the area containing regression trees and is expressed as: (4) $F=\left\{f\left(x\right)={\omega }_{q\left(x\right)}\right\}\left(q:R^n\to T,\omega \in R^T\right)$(4) where q indicates the associated leaf for each data point as a tree structure, T is the number of leaves on the tree, f_p is a function corresponding to a stand-alone tree structure, and ω is the output weight of leaves.

The target function (J) for the XGBoost method is defined as in Eq. (5) to minimize the error of ensemble trees: (5) $J^{\left(t\right)}=\sum _{i=1}^{m}L\left(y_i,{\hat{y}}_i^{\left(t-1\right)}+f_t\left(x_i\right)\right)+\mathrm{\Omega }\left(f_t\right)$(5) where L is the training loss function used to determine the distance between the estimated (${\hat{y}}_i$) and observed (y_i) values, superscript t denotes the number of iterations, and Ω is the regularization function that controls the model complexity: (6) $\mathrm{\Omega }\left(f\right)=\gamma T+\frac{1}{2}\lambda {\omega }^2$(6) where λ is the penalty coefficient, and γ shows the complexity of each tree. A higher value of γ means fewer complex trees. The hyperparameters of XGBoost are the learning rate, max_delta_step, max_depth, and min_child_weight. The learning rate demonstrates the shrinkage (reduction in the size of incremental steps) used to prevent overfitting in the learning process. A shrinking weight in each step will result in a more conservative model. Max_delta_step indicates the weight estimate of each tree. A positive value results in a conservative update step. The maximum depth of the tree, max_depth, is used to control overfitting in the learning phase. Min_child_weight is the minimum summation of weights in a child tree, higher values of which make the model more conservative. In this study, we apply three optimization techniques to find the optimal hyperparameters of XGBoost, as described in the following subsections.

3.2. PSO

PSO, introduced by Kennedy and Eberhart (1995), is a metaheuristic algorithm that draws inspiration from nature. Each particle has an assigned position (x) and velocity (v). The initial position (x) and velocity (v) of each particle are defined by an n-dimensional random vector and vector of zeros, respectively. Each dimension represents a specific decision parameter. Assume that particle i moves at velocity $v_i^t$ at iteration t, as a part of a swarm sized P, where $i\in \left\{1,2,\dots ,P\right\}$. The particle velocity at the time $t+1$ ($v_i^{t+1}$.) can be obtained as (Le et al., 2019): (7) $v_i^{t+1}=w{v}_i^t+c_1{r}_1\left(x_i^{\ast }-x_i^t\right)+c_2{r}_2\left(g-x_i^t\right)$(7) where $x_i^{\ast }$ represents the best location achieved by particle i upo the current iteration, g is the present best global solution among all swarms w is the inertial weight, c₁ and c₂ are the cognitive and social parameters, respectively; and r₁ and r₂ are random numbers that vary from 0 to 1.

In the following step, given the updated velocity of particle i, the particle's position is adjusted as follows: (8) $x_i^{t+1}=x_i^t+v_i^{t+1}$(8) where $x_i^t$indicates the location of particle i at iteration t.

3.3. RFO

Red foxes are highly populated species of foxes that can live and survive in different climatic conditions. Based on the red foxes’ lifestyle and method of hunting, Połap and Woźniak (2021) developed a metaheuristic optimization algorithm named RFO. The RFO initialization is simulated by a fixed number of foxes, with the coordinates of each fox defined as: (9) $X=\left[x_0,x_1,\dots ,x_{n-1}\right]$(9) where n describes the number of coordinates.

Each fox in each iteration is identified using the notation $(X_j^i{)}^t$, where i indicates the number of a fox in the population, j is the coordinate based on the dimension of the searching area, and t is the number of iterations.

Assuming $f\in R^n$ to be a condition function of n variables, each point in the search space $[a,b{]}^n$, where $a,b\in R$, can be expressed as in Eq. (10): (10) $(X{)}^i=[(x_0{)}^i,(x_1{)}^i,\dots ,(x_{n-1}{)}^i]$(10) (X)ⁱ is considered the optimum solution once the value of f ((X)ⁱ) yields the global optimal. Every fox must participate in protecting the flock from threats. If adequate prey is not available in the area, the foxes must move to farther regions to achieve a better result for the exploration term. The information collected by a fox for the best location is shared with the others if a region with adequate prey is found. The fitness value is used to sort the results of the evaluated function.

The Euclidean distance (D) is used to find $(X^{best}{)}^t$: (11) $D\left[{\left(X^i\right)}^t,{\left(X^{best}\right)}^t\right]=\sqrt{||{\left(X^i\right)}^t-{\left(X^{best}\right)}^t||}$(11) where || . || is the Euclidean norm. Thus, the foxes travel toward the optimum solution as follows: (12) $(X^i{)}^t=(X^i{)}^t+a\times sign((X^{best}{)}^t-(\left(X{)}^i{)}^t\right)$(12) where a is a randomly chosen integer number in the range of $a\in \left(0,D\left[{\left(X^i\right)}^t,{\left(X^{best}\right)}^t\right]\right)$

Then, the appropriate solution is obtained using the updated position of the candidates. Otherwise, the former location is retained. The probability of a fox being observed while they move to capture the prey is modeled using a random variable λ: (13) $\{\begin{array}{l}movecloser\lambda >0.75\\ stayandcamouflage\lambda \le 0.75\end{array}$(13) The radius (r) is another important term in the RFO, which indicates the vision radius of the predator fox: (14) $r=\{\begin{array}{l}\alpha {\displaystyle \frac{\sin {\phi }_0}{{\phi }_0}}{\phi }_0\ne 0\\ \text{β}{\phi }_0=0\end{array}$(14) where α is the scaling coefficient for modeling changes in the fox’s vision radius while they approach a prey and lies in the range [0, 0.2]. β is a random value between 0 and 1 that denotes the impact of adverse weather such as rain or fog on the foxes’ vision angle. φ₀, which represents the fox’s angle of sight, varies in the range [0, 2π]. Therefore, the fox population approaching prey can be simulated as follows: (15) $\begin{array}{rl}& \{\begin{array}{l}{x}_0^{New}=\alpha r\cos \left({\phi }_1\right)+x_0^{Old}\\ x_1^{New}=\alpha r\sin \left({\phi }_1\right)+\alpha r\cos \left({\phi }_2\right)+x_1^{Old}\\ x_2^{New}=\alpha r\sin \left({\phi }_1\right)+\alpha r\sin \left({\phi }_2\right)+\alpha r\cos \left({\phi }_3\right)+x_2^{Old}\\ ...\\ x_{n-2}^{New}=\alpha r\sum _{k=1}^{n-2}\sin \left({\phi }_k\right)+\alpha r\cos \left({\phi }_{n-1}\right)+x_{n-2}^{Old}\\ x_{n-1}^{New}=\alpha r\sin \left({\phi }_1\right)+\alpha r\sin \left({\phi }_2\right)+...+\alpha r\sin \left({\phi }_{n-1}\right)\\ +x_{n-1}^{Old}\end{array}\end{array}$(15) where x^New is the updated coordinate, x^Old is the coordinate in the previous iteration, and ${\phi }_1,{\phi }_2,\dots .{\phi }_{n-1}$ are the foxes’ angles of sight for the corresponding coordinates.

3.4. RPSO

RPSO is a variant of the PSO, which was introduced by Roder et al. (2020). We introduce the theory of relativity to illustrate the RPSO. The theory of relativity is one of the most significant theories in physics, proposed by Albert Einstein in 1916. Conventional mechanics’ theories were implemented to numerically model the phenomena detected by Einstein. One of the most famous theories pertains to momentum, which can be calculated in three-dimensional coordinates as: (16) $M\left(v\right)=\varsigma \left(v\right)mv$(16) where M is the momentum, v = (v_x, v_y, v_z) and m are the velocity and mass of a body, respectively. $\varsigma$ is the Lorentz factor, which is defined as follows: (17) $\varsigma \left(v\right)=\frac{1}{\sqrt{1-{\left(\frac{|v|}{c}\right)}^2}}$(17) where |v| denotes the size of vector v, and c is the speed of light.

Unlike PSO, RPSO considers the effects of the speed of light, improved social behaviour of swarms, and particles’ mass on the optimum values of the solution. RPSO transforms the three-dimensional momentum formula (Eq. 16) into an n-dimensional formula to compute the velocity of each particle in an n-dimensional search space: (18) $v_i^{t+1}=M\left(v_i^t\right)+c_1{r}_1\left(x_i^{\ast }-x_i^t\right)+c_2{r}_2\left(g-x_i^t\right)$(18) where M(.) is the typical relativistic momentum in three dimensions, which is extended to the n-dimensional space. Mass m is selected from a uniform distribution with values between [0, 1] to compute the relativistic momentum of each particle (Roder et al., 2020).

3.5. Model performance evaluation metrics

Three statistical indices are used to assess the efficiency of the proposed methods: root mean square error (RMSE), mean absolute error (MAE), and coefficient of determination (R²). Lower MAE and RMSE denote better results. R², which indicates the percentage of variation in dependent variable that can be accounted for by independent variables, lies in the range [0,1]. R² values closer to 1 indicate a stronger relationship between the predictions and observations. The three metrics are defined as follows: (19) $\begin{array}{rl}\text{RMSE}& =\sqrt{\frac{{\sum }_{i=1}^{\text{N}}{\left(d_{s{e}_i}^{obs}-d_{s{e}_i}^{pred}\right)}^2}{\text{N}}}\end{array}$(19) (20) $\begin{array}{rl}\text{MAE}& =\frac{{\sum }_{i=1}^{\text{N}}|d_{s{e}_i}^{obs}-d_{s{e}_i}^{pred}|}{\text{N}}\end{array}$(20) (21) $\begin{array}{rl}{\text{R}}^2& =1-\frac{{\sum }_{i=1}^{\text{N}}{\left(d_{s{e}_i}^{obs}-d_{s{e}_i}^{pred}\right)}^2}{{\sum }_{i=1}^{\text{N}}{\left(d_{s{e}_i}^{obs}-{\overline{d}}_{se}\right)}^2}\end{array}$(21) where N is the number of data points, $d_{s{e}_i}^{pred}$ and $d_{s{e}_i}^{obs}$ are the estimated and observed scour depths for the ith data point, respectively, and ${\overline{d}}_{se}$ indicates the average values of observed scour depth.

3.6. SHAP

Lundberg and Lee (2017) developed the SHAP method based on game theory to interpret ML models. The concept was used in game theory to assess a player’s participation in a cooperative team game and distribute a fair reward among players based on their contributions (Shapley, 1953). SHAP determines the degree of contribution of each input to the model’s output prediction. Moreover, it evaluates whether the effect of each input on the model’s prediction is positive or negative. The use of SHAP to interpret ML models can help prevent the problem of these models operating as black boxes. The model behaviour in both global (for the entire dataset) and local (for a single prediction) aspects is explained by assigning SHAP values to each feature (Zhang et al., 2022). The contribution of each input feature can be determined as: (22) ${\kappa }_i=\sum _{\zeta \subseteq G}\frac{|\zeta |!\left(n-|\zeta |-1\right)!}{n!}\left[\tau \left(\zeta \cup \left\{i\right\}\right)-\tau \left(\zeta \right)\right]$(22) where ${\kappa }_i$ indicates the contribution of the i-th feature (in terms of Shapley values), G is the collection of all features, n is the number of features in G, and $\zeta$ is a subset of G, |.| denotes the number of members in a set, $\tau \left(\zeta \cup \left\{i\right\}\right)$ and $\tau \left(\zeta \right)$indicate the mel’s output with and without the i-th feature, respectively. The model output can be written as a linear summation of a constant base value and SHAP values as (Khattak et al., 2022): (23) $e\left(z^{\mathrm{\prime }}\right)={\kappa }_0+\sum _{i=1}^n{\kappa }_i{z}^{\mathrm{\prime }}$(23) where ${\kappa }_0$ is a base value, $z^{\mathrm{\prime }}\in \{0,1{\}}^n$, where $z^{\mathrm{\prime }}=1$ when a feature is present and $z^{\mathrm{\prime }}=0$otherwise.

Figure 1 illustrates the process of developing hybrid ML models for scour depth prediction.

Figure 1. Process flow of the proposed hybrid ML models and SHAP for estimating the scour depth near circular bridge piers.

4. Results and discussion

4.1. Model evaluation

The 841 experimental data points are split randomly into training (70%), validation (15%), and testing (15%) datasets. The training dataset is used to train the ML models and typically includes most of the data. The model uses the training data to minimize the differences between observed and predicted values by adjusting the hyperparameters (Brownlee, 2020). The validation dataset is used to assess the model efficiency during training and to obtain the best set of hyperparameters. The validation dataset is distinct from the training dataset to prevent overfitting (Brownlee, 2020). The performance of the final (calibrated) ML model is assessed using the testing dataset, which is completely withheld from the model during the training and validation stages (Mangalathu et al., 2020). In this manner, testing data enable an unbiased evaluation of ML models. All three datasets should be used to ensure that the ML model is precise, robust, and not overfitted to training data.

In this study, the hyperparameters of the XGBoost method, i.e. learning rate, max_delta_step, max_depth, min_child weight, and n_estimators are optimized using PSO, RFO, and RPSO. The default values of these parameters are used to initialize XGBoost. The merits of hybrid models include higher accuracy, flexibility, and robustness in handling high-dimensional data (Jalil et al., 2022). Tables 3 and 4 present the optimal values of the hyperparameters obtained by different optimization approaches for dimensional and dimensionless datasets, respectively.

Table 3. Initial hyperparameters of the XGBoost model and values optimized using the PSO, RFO, and RPSO techniques for the dimensional dataset.

Parameter	XGBoost	PSO–XGBoost	RFO–XGBoost	RPSO–XGBoost
learning_rate	0.3	0.7917	0.1350	0.2575
max_delta_step	0	9.7905	1.0746	9.9969
max_depth	6	10	19	17
min_child_weight	1	1.1900	6.0487	1.1051
Number of iterations	100	72	369	369

Table 4. Initial hyperparameters of the XGBoost model and values optimized using the PSO, RFO, and RPSO techniques for the dimensionless dataset.

Parameter	XGBoost	PSO–XGBoost	RFO–XGBoost	RPSO–XGBoost
learning_rate	0.3	0.4902	0.3994	0.4848
max_delta_step	0	8.0483	7.1944	2.9955
max_depth	6	17	20	6
min_child_weight	1	7.142	10	7.69
Number of iterations	100	181	188	363

First, the original variables in Eq. (1) are used in the developed models to estimate the scour depth. Figure 2 presents the scatter plot of scour depth predictions from XGBoost, PSO–XGBoost, RFO–XGBoost, and RPSO–XGBoost versus observations for the training, validation, and testing stages. The scour depth estimations from all models are close to the 45° line, which shows that the proposed models can efficiently predict the scour depth. However, the predictions for larger scour depths deviate from the 1:1 line in the validation and testing steps, owing to the lower frequency of data in this range. Scour predictions in the training and validation phases are better than those in the testing step because the training and validation data are used in the learning process.

Figure 2. Scatter plot of predicted versus observed scour depth values for the training, validation, and testing phases using XGBoost, PSO–XGBoost, RFO–XGBoost, and RPSO–XGBoost.

Table 5 summarizes the performance metrics of the proposed models (i.e. RMSE, MAE, and R²) for the training, validation, and testing stages.

Table 5. Statistical indicators of scour depth predictions made by the proposed machine learning techniques using the dimensional variables in Eq. (1).

	Training data			Validation data			Testing data
Method	RMSE (m)	MAE (m)	R^2	RMSE (m)	MAE (m)	R^2	RMSE (m)	MAE (m)	R^2
XGBoost	0.0102	0.0035	0.9935	0.0255	0.0154	0.8644	0.0434	0.0228	0.9152
PSO–XGBoost	0.0101	0.0033	0.9936	0.0191	0.0132	0.9238	0.0424	0.0239	0.9188
RFO–XGBoost	0.0101	0.0029	0.9936	0.0171	0.0113	0.9392	0.0423	0.0209	0.9191
RPSO–XGBoost	0.0101	0.0023	0.9936	0.0151	0.0104	0.9526	0.0398	0.0207	0.9284

As expected, XGBoost is the least accurate method with MAEs of 0.0035, 0.0154, and 0.0228 m and RMSEs of 0.0102, 0.0255, and 0.0434 m in the training, validation, and training stages, respectively. Figure 2 and Table 5 demonstrate that all three hybrid models are more accurate than the XGBoost method owing to the use of optimized hyperparameters. In the training phase, PSO–XGBoost, RFO–XGBoost, and RPSO–XGBoost reduce the RMSE of XGBoost by 0.98%. The corresponding RMSE reductions are 2.30%, 2.53%, and 8.29% for testing and 25.09%, 32.94%, and 40.78% for validation. In terms of the MAE, in the training step, PSO–XGBoost, RFO–XGBoost, and RPSO–XGBoost enhance the accuracy of XGBoost by 5.71%, 17.14%, and 34.28%, respectively. RPSO–XGBoost exhibits the highest performance with the lowest RMSEs (0.0101, 0.0151, and 0.0398 m) and MAEs (0.0023, 0.0104, and 0.0207 m) for the training, validation, and testing phases, respectively, followed by RFO–XGBoost and PSO–XGBoost. In the training step, the MAEs of RPSO–XGBoost are 30.30% and 20.68% lower than those of PSO–XGBoost and RFO–XGBoost, respectively. In the validation step, RPSO–XGBoost improves the MAE (RMSE) of PSO–XGBoost and RFO–XGBoost by 21.21% (40.78%) and 7.96% (11.69%), respectively. A similar trend of MAE (RMSE) reduction is seen in the testing stage, where the MAEs (RMSEs) of RPSO–XGBoost are 13.38% (6.13%) and 0.95% (5.91%) smaller than those of PSO–XGBoost and RFO–XGBoost, respectively. Overall, the proposed models can be ranked in the decreasing order of their ability to predict the scour depth around bridge piers as RPSO–XGBoost, RFO–XGBoost, PSO–XGBoost, and XGBoost.

Subsequently, the dimensionless parameters in Eq. (2) are used in the models to estimate d_se/Y. Figure 3 compares the d_se/Y estimates from XGBoost, PSO–XGBoost, RFO–XGBoost, and RPSO–XGBoost with observations for the training, validation, and testing stages.

Figure 3. Scatter plot of predicted versus observed dimensionless scour depth (d_se/Y) in the training, validation, and testing phases, derived using XGBoost, PSO–XGBoost, RFO–XGBoost, and RPSO–XGBoost.

Figure 3 shows that the scour depth predictions from all proposed models are concentrated around the 45° line, indicating their satisfactory accuracies. A comparison of Figs. 2 and 3 demonstrates that all models produce more reliable results when trained by dimensional data (Eq. (1)) instead of dimensionless data (Eq. (2)). A possible explanation for this may be that the use of raw variables rather than a combination of variables enhances the flexibility of models in simulating the highly nonlinear relationship between the inputs and output. This result is in line with those of Bateni et al. (2007a, 2007b) and Zounemat-Kermani et al. (2009), who reported that the results of ANN and ANFIS models trained by the dimensional dataset were superior to the models trained by dimensionless data.

Table 6 presents a comparison of the predictive performances of the proposed methods for d_se/Y. The MAEs (RMSEs) of the hybrid models are lower than those of XGBoost. The best d_se/Y estimates are attributed to RPSO–XGBoost with RMSE = 0.1723, MAE = 0.0431, and R² = 0.9765 for training; RMSE = 0.2065, MAE = 0.1354, and R²= 0.9714 for validation; and RMSE = 0.2065, MAE = 0.1354, and R²= 0.9714 for testing. For the testing phase, the RMSE of RPSO–XGBoost is 4.36% and 5.48% lower than those of RFO–XGBoost and PSO–XGBoost, respectively. XGBoost has the weakest performance among all proposed methods with RMSE = 0.1825, MAE = 0.0763, and R² = 0.9735 in the training phase; RMSE = 0.4066, MAE = 0.1832, and R² = 0.8889 in the validation phase; and RMSE = 0.3435, MAE = 0.1712, and R² = 0.7824 in the testing phase. Table 6 and Figure 3 also highlight the superiority of all hybrid models over XGBoost: In the training stage, the MAEs (RMSEs) of PSO–XGBoost, RFO–XGBoost, and RPSO–XGBoost are 14.94%, (0.11%), 30.27% (3.29%), and 43.51% (5.59%) lower than that of XGBoost, respectively. Similarly, for validation data, the MAEs (RMSEs) of PSO–XGBoost, RFO–XGBoost, and RPSO–XGBoost are lower than those of XGBoost by 5.40% (20.76%), 19.05% (39.69%), and 26.09% (49.21%), respectively. These results indicate that PSO, RFO, and RPSO can successfully tune the hyperparameters of XGBoost, thereby enhancing its scour depth prediction capabilities. Therefore, the proposed models can be ranked in decreasing order of their accuracy in predicting d_se/Y as RPSO–XGBoost, RFO–XGBoost, PSO–XGBoost, and XGBoost. This finding is consistent with the previous ranking of methods for scour depth prediction using the dimensional variables in Eq. (1)(1) $d_{se}=f\left(\rho ,\mu ,V,Y,g,d_{50},V_c,D\right)$(1) .

Table 6. Statistical indices of (d_se/Y) predictions from the proposed methods using the dimensionless variables in Eq. (2).

	Training data			Validation data			Testing data
Method	RMSE	MAE	R^2	RMSE	MAE	R^2	RMSE	MAE	R^2
XGBoost	0.1825	0.0763	0.9735	0.4066	0.1832	0.8889	0.3435	0.1712	0.7824
PSO–XGBoost	0.1823	0.0649	0.9736	0.3222	0.1733	0.9302	0.3340	0.1690	0.7945
RFO–XGBoost	0.1765	0.0532	0.9753	0.2452	0.1483	0.9596	0.3301	0.1593	0.7992
RPSO–XGBoost	0.1723	0.0431	0.9765	0.2065	0.1354	0.9714	0.3157	0.1657	0.8164

The superiority of RPSO–XGBoost over PSO–XGBoost and RFO–XGBoost is attributed to several factors. Firstly, the RPSO technique incorporates relativistic effects, allowing particles to move at speeds that can approach the speed of light (Roder et al., 2020). These relativistic effects enable RPSO–XGBoost to explore the search space extensively, conducting a more comprehensive and thorough search for optimal hyperparameters (Roder et al., 2020). Secondly, the relativistic effects in RPSO facilitate more efficient convergence of particles towards the global optimum (Roder et al., 2020). The effective navigation of the complex search space by RPSO allows RPSO–XGBoost to reveal significant patterns and dependencies that influence scour depth, resulting in highly accurate predictions.

Based on testing data, Table 7 compares the performances of the proposed models and 28 existing techniques. The numbers in parentheses show the percentage of increase in the RMSE/MAE and decrease in R² compared with those of RPSO–XGBoost, which is the highest-performing approach developed in this study. Because the testing data are withheld from the model in the training and validation steps, they can be used as an independent data source to perform a fair evaluation of the final models. The results show that the four proposed methods outperform all the existing approaches. Among the 28 existing methods, the model proposed by Shamshirband et al. (2020) achieves the best scour depth predictions with RMSE, MAE, and R² of 0.0654, 0.0470, and 0.8268 m, respectively.

Table 7. Comparison of scour depth estimates from XGBoost, PSO–XGBoost, RFO–XGBoost, and RPSO–XGBoost with those obtained using existing approaches.

Methods	RMSE (m)	MAE (m)	R^2
Inglis-Poona I (Inglis, 1949)	0.1272 (68.71%)	0.0892 (76.79%)	0.3617 (156%)
Inglis-Poona II (Inglis, 1949)	0.0921 (56.78%)	0.0609 (66%)	0.7872 (17.93%)
Laursen and Toch (1956)	0.1009 (60.55%)	0.0580 (64.31%)	0.7522 (23.42%)
Chitale (1962)	0.2090 (80.95%)	0.1220 (83.03%)	0.1002 (826%)
Blench-Inglis I (Blench et al., 1962)	0.0885 (55.03%)	0.0604 (65.72%)	0.7830 (18.56%)
Blench-Inglis II (Blench et al., 1962)	0.1423 (72.03%)	0.0791 (73.83%)	0.2343 (296%)
Larras (1963)	0.1000 (60.2%)	0.0668 (69.01%)	0.7158 (29.70%)
Laursern I (Neill, 1964)	0.0926 (57.01%)	0.0634 (67.35%)	0.8727 (6.38%)
Breusers (1965)	0.1018 (60.90%)	0.0552 (62.5%)	0.7334 (26.58%)
Shen et al. (1969)	0.1007 (60.47%)	0.0710 (70.84%)	0.7471 (24.26%)
Blench (1969)	0.7712 (94.83%)	0.7463 (97.22%)	0.1318 (604%)
Shen (1971)	0.0828 (51.93%)	0.0586 (64.67%)	0.7321(26.81%)
Hancu (1971)	0.1051 (62.13%)	0.0714 (71%)	0.5507 (68.58%)
Coleman (1971)	0.1889 (78.93%)	0.1309 (84.18%)	0.7368 (26%)
Neill (1973)	0.1137 (64.99%)	0.0595 (65.21%)	0.7334 (26.58%)
Breusers et al. (1977)	0.1053 (62.20%)	0.0590 (64.91%)	0.5539 (67.61%)
Jain and Fischer (1979)	0.0784 (49.23%)	0.0566 (63.42%)	0.8561 (8.44%)
Froehlich (1988)	0.1360 (70.73%)	0.0865 (76.07%)	0.8488 (9.37%)
Froehlich design equation (1988)	0.0924 (56.92%)	0.0560 (63.03%)	0.7832 (18.54%)
Melville and Sutherland (1988)	0.1388 (71.32%)	0.0854 (75.76%)	0.8706 (6.64%)
Ansari and Qadar (1994)	0.1827 (78.21%)	0.1307 (84.16%)	0.6828 (35.96%)
Simplified Chinese equation (Landers & Mueller, 1996)	0.1283 (68.97%)	0.0678 (69.46%)	0.3873 (139%)
Melville (1997)	0.0984 (59.55%)	0.0566 (63.42%)	0.7633 (21.63%)
Lyn (2008)	0.1739 (77.11%)	0.1090 (81%)	0.1095 (747%)
HEC-18 (Arneson et al., 2012)	0.0763 (47.83%)	0.0552 (62.5%)	0.8078 (14.93%)
Sheppard et al. (2014)	0.0742 (46.36%)	0.0471 (56.05%)	0.8443 (9.96%)
Choi et al. (2017)	0.2552 (84.40%)	0.1890 (89.04%)	0.4278 (117%)
Shamshirband et al. (2020)	0.0654 (39.14%)	0.0470 (55.95%)	0.8268 (12.28%)
XGBoost (this study)	0.0434	0.0228	0.9152
PSO–XGBoost (this study)	0.0424	0.0239	0.9188
RFO–XGBoost (this study)	0.0423	0.0209	0.9191
RPSO–XGBoost (this study)	0.0398	0.0207	0.9284

Shamshirband et al. (2020) used the PSO algorithm to develop an explicit equation for estimating scour depth around circular bridge piers. Their equation was developed by minimizing the discrepancy between predicted and observed scour depth values via the PSO algorithm (Shamshirband et al., 2020). The PSO technique enhances the capability of their expression to capture inherent relationships among input and output data, leading to improved prediction accuracy. Furthermore, the optimization process performed by the PSO algorithm facilitates the generalizability of the equation based on the available dataset (Taormina & Chau, 2015). As a result, their derived equation has the potential to be applied to a diverse range of scenarios beyond the specific dataset utilized for its development. The ability to generalize effectively is of considerable significance for practical applications.

Other studies in Table 7 mostly used regression techniques to derive equations for predicting d_se. The traditional regression-based approaches cannot robustly capture the nonlinearity and complexity among the scour depth and its influential variables (Ebtehaj et al., 2017).

The RMSE (MAE) of the scour depth estimates from RPSO–XGBoost is 39.14% (55.95%) lower than those of Shamshirband et al. (2020). The model proposed by Blench (1969) produces the least accurate results, with RMSE = 0.7712 m, MAE = 0.7463 m, and R² = 0.1318. On average, RPSO–XGBoost improves the RMSE and MAE of scour depth predictions by 64.91% and 70.94% compared with those of the 28 existing approaches. Figure 4 visually compares the scour depth estimates from the four proposed models with those obtained using the model of Shamshirband et al. (2020) as the highest-performing method among the 28 considered techniques. Most of the scour depth predictions obtained using the proposed models are closer to the 1:1 line, which demonstrates their superiority over those of the model proposed by Shamshirband et al. (2020).

Figure 4. Visual comparison of the predictions obtained using the proposed methods and the method proposed by Shamshirband et al. (2020), based on testing data

The violin plot in Fig. 5a shows the distribution of scour depth estimates from RPSO–XGBoost (as the best proposed model) in comparison with those obtained using two famous models (i.e. HEC-18 and Sheppard et al. (2014)) and the best method among the 28 considered existing methods (i.e. Shamshirband et al. (2020)). A similar comparison is made for d_se/Y in Fig. 5b. In general, violin plots combine a kernel density plot with a box plot. Unlike box plots, violin plots can show both the statistic and density of data, thereby explaining the variability of data. The small white dot in the violin plots in Fig. 5 represents the median of the data. The interquartile range is presented by a thick black line. The thin black line indicates the data beyond the interquartile range except outliers. The shape of the data distribution is shown by the kernel density approximation on the two sides of the black line. The wider shape of the violin plot around the median denotes a high concentration of data in this region. The tapered shape of the ends of the violin plot indicates a lower concentration of data in that area.

Figure 5 shows that the median and interquartile values of the scour depths predicted by the RPSO–XGBoost are close to the observations, highlighting its high predictive accuracy. For both dimensional and dimensionless data, the scour depth estimates from all methods are clustered around the median, indicating a higher probability in this region. In the case of dimensional data, the distribution of scour depth predictions from RPSO–XGBoost is similar to that of the measurements. In addition, the interquartile range, range of data, and median of predictions from RPSO–XGBoost are closer to the observations for both dimensional and dimensionless data, compared with those of Shamshirband et al. (2020), Sheppard et al. (2014), and HEC-18. A comparison of Figure 5(a,b) confirms the findings of the previous analyses demonstrating the superiority of scour depth predictions obtained using the dimensional variables in Eq. (1) over those obtained using the dimensionless variables in Eq. (2). The upper end of the violin plots for d_se/Y predictions (Figure 5(b)) indicates a higher value than that of d_se/Y measurements because the models overestimate d_se/Y in the case of large values.

Figure 5. Violin plots comparing the distribution of (top) dimensional dataset and (bottom) dimensionless dataset estimates from different approaches with measurements.

To further investigate the ability of the proposed methods, the Taylor diagram for scour depth measurements and estimates is derived, as shown in Figure 6. Taylor diagrams are valuable tools for visualizing results, as they can graphically summarize several evaluation metrics, i.e. the SD, correlation coefficient between the predicted and measured values, and centered RMSE, in a single plot. Thus, Taylor diagrams can help obtain reliable conclusions regarding the efficiency of the proposed models. Notably, the statistical indicators shown in the Taylor diagram in Figure 6 correspond to the testing data. When using the dimensional dataset, the proposed XGBoost, PSO–XGBoost, and RPSO–XGBoost models underestimate the range of scour depth variation, whereas RFO–XGBoost tends to overestimate the scour depth variations. When using dimensionless data, all the proposed methods overestimate the variations in the scour depth, as indicated by the higher SD compared with the observations (i.e. the values lie beyond the red dashed curve). Additionally, Figure 6 provides further evidence of the higher accuracy of proposed models using dimensional variables compared with those using dimensionless variables, as the correlation coefficient of the former is more than 0.95 but that of the latter ranges between 0.9 and 0.95. RPSO–XGBoost yields better results for both input configurations given their higher correlation with the observations and lower RMSE compared with those of the other three methods. Furthermore, the variability (represented by SD) of d_se and d_se /Y predicted by all proposed methods is close to that of the observations. For both dimensional and dimensionless datasets, the predictions of XGBoost, PSO–XGBoost, RFO–XGBoost, and RPSO–XGBoost are closer to the scour depth observations (black triangle in the Taylor diagram) compared with those obtained using the method of Shamshirband et al. (2020), Sheppard et al. (2014), and HEC-18 (Arneson et al. 2012). Furthermore, the SD of d_se/Y predictions from Sheppard et al. (2014) and HEC-18 differs significantly from the observations, which shows that these methods cannot accurately predict the variability in d_se/Y.

Figure 6. Taylor diagrams for scour depth estimates from various approaches and observations.

Figure 7(a,b) show the RMSEs of d_se and d_se/Y predictions from RPSO–XGBoost (the most accurate method) under different numbers of iterations, respectively. As the number of iterations increases, the RMSEs of d_se and d_se/Y predictions from RPSO–XGBoost decrease rapidly to their asymptotes. A learning curve of a model is considered a good fit (i.e. neither underfit nor overfit) if it satisfies two conditions: (1) The RMSEs in the training phase decrease to the final stable value, and (2) the RMSEs in the validation phase follows the trend of the training phase with minor differences. According to Figure 7, the learning behaviour of the proposed model (RPSO–XGBoost) can be a good fit.

Figure 7. Variations in the RMSE of (left) d_se and (right) d_se/Y predictions from RPSO–XGBoost (the best model in this study) versus the number of iterations in the training and validation stages.

4.2. Model interpretation

A SHAP analysis is performed to identify the variable with the most notable contribution to the scour depth predictions. As mentioned previously, SHAP values not only allow us to determine the contribution of each variable contributes to the prediction but also to identify and visualize the significant relationships in the model. SHAP values can be visualized using various approaches such as waterfall, force plot, decision plot, mean SHAP plot, and beeswarm plot. A beeswarm plot aggregates the SHAP values of all observations to ensure that all SHAP values can be visualized at once. Figure 8(a,b) show the beeswarm plots based on RPSO–XGBoost (the most accurate model) for dimensional and dimensionless variables, respectively.

Figure 8. SHAP summary plots based on RPSO–XGBoost to rank the significance of each feature in estimating (a) d_se and (b) (d_se/Y).

In this plot, the variables are ordered based on their importance on the y-axis. The x-axis denotes the mean SHAP value, and the colour of each point shows the actual value of that feature (i.e. not the SHAP value). Thicker parts of the beeswarm plot imply a higher density of predictions in that region. Figure 8(a) indicates that for dimensional data, the pier diameter (D) has the most notable effect on scour depth prediction. The second most influential variable is the flow depth (Y), followed by the flow velocity (V), median grain size (d₅₀), and critical velocity of sediments (V_c). This ranking is consistent with the study of Bateni et al. (2007a), who reported that the pier diameter and critical velocity have the greatest and least impact on the estimation of scour depth using ANN, respectively. Larger values of the pier diameter correspond to larger SHAP values, which in turn yield higher scour depth estimates. A similar positive relationship is observed for the flow depth and flow velocity. In contrast, a clear inverse relationship is found between the scour depth and the median grain size and critical velocity, both of which are associated with sediment characteristics. These relationships are compatible with the physical nature of the scouring process, as larger pier diameter, flow depth, and flow velocity result in higher scour depths, and piers located in riverbeds with finer sediments and sediments with a lower critical velocity typically experience larger scour depths. As outlined in Figure 8(b), the most and least significant dimensionless variables for predicting d_se/Y are D/Y and Re, respectively. The remaining dimensionless variables can be listed in decreasing order of their importance as V/V_c, d₅₀/Y, and Fr. Higher values of D/Y and V/V_c lead to higher d_se/Y predictions. Because of certain outliers in the SHAP values of d₅₀/Y, Fr, and Re, no clear relationship can be concluded between these variables and d_se/Y predictions.

To ensure the physical consistency of scour depth estimates from ML models, the following multifaceted approach is essential (Najafzadeh & Oliveto, 2020): (1) Compare predictions from ML models with those obtained from established physical models and/or widely used empirical equations to verify the physical consistency of results (Kollet et al., 2017). In this study, we met this criterion by comparing the scour depth estimations from XGBoost, PSO–XGBoost, RFO–XGBoost, and RPSO–XGBoost with those of 28 existing equations in the literature (see Table 7); (2) Identify the importance of input features to ensure an alignment with the physical governing principles. Herein, the SHAP method is used to find the relative significance of each input feature on scour depth predictions. The results of SHAP are consistent with those of Bateni et al. (2007a, 2007b), Najafzadeh and Barani (2011), Sharafi et al. (2016), Ebtehaj et al. (2019), and Shamshirband et al. (2020); (3) Evaluate the performance of ML models using unseen data to provide additional evidence of their physical consistency and generalization capability. This study assessed the feasibility of all ML models via the testing dataset. This dataset is kept hidden from the models during the training and validation phases. Hence, this study considered all the criteria to ensure the physical consistency of our results.

5. Conclusions

Three novel hybrid ML methods: PSO–XGBoost, RFO–XGBoost, and RPSO –XGBoost are applied to estimate the equilibrium scour depth around circular bridge piers. A comprehensive laboratory and field dataset containing 841 data points from 35 studies is used to train, validate, and test the proposed ML methods. The hyperparameters of the XGBoost method are tuned using three optimization techniques (PSO, RFO, and RPSO). The proposed models are trained using both dimensional and dimensionless datasets. The results are evaluated using three statistical metrics i.e. the RMSE, MAE, and R². All three proposed methods accurately predict the scour depth. The higher accuracy of the hybrid models over XGBoost shows that the three optimization techniques can efficiently tune the hyperparameters of XGBoost and enhance its accuracy in scour depth prediction.

For models trained by dimensional data, RPSO–XGBoost yields the best scour depth predictions with RMSE = 0.0101 m, MAE = 0.0023 m, and R² = 0.9936 in training; RMSE = 0.0151 m, MAE = 0.0104 m, and R² = 0.9526 in validation; and RMSE = 0.0398 m, MAE = 0.0207 m, and R² = 0.9284 in the testing phrase. Similarly, for dimensionless data, RPSO–XGBoost produces the most accurate scour depth predictions with RMSE = 0.1723, MAE = 0.0431, and R² = 0.9765 in training; RMSE = 0.2065, MAE = 0.1354, and R² = 0.9714 in validation; and RMSE = 0.3157, MAE = 0.1657, and R² = 0.8164 in testing. Moreover, more accurate results can be achieved when the proposed models are trained with dimensional data rather than dimensionless data.

The results of the proposed methods are compared with those of 28 existing methods, and their superiority over the existing techniques is demonstrated. The MAE (RMSE) of RPSO–XGBoost (as the best method from our study) is 55.95% (39.14%) lower than that of the model proposed by Shamshirband et al. (2020) (as the best method among the 28 studies). Finally, the aggregated (SHAP) values of each variable for both dimensional and dimensionless datasets are used to determine the parameters, most notably affecting the scour depth predictions and clarify the input–output relationships. The SHAP analysis shows that among the dimensional variables, the pier diameter (D) and critical velocity (V_c) have the most and least notable effects on the scour depth estimation, respectively. Higher values of D, flow depth (Y), and flow velocity (V) yield larger values of the scour depth, whereas higher median grain size (d₅₀) and V_c lead to lower values of the scour depth. SHAP analysis on the dimensionless data shows that D/Y and Re have the most and least notable contributions to d_se/Y, respectively.

This finding highlights the marginal impact of Re on the scour depth around circular bridge piers.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This study has been made possible by the Hawaii Department of Transportation (HDOT) and Federal Highway Administration (FHWA) grant DOT-10-030 to the University of Hawaii at Manoa.

Appendix

Table A presents several well-known equations to determine the depth of scour near bridge piers. In the equations, d_se is the equilibrium scour depth, Y is the flow depth, D is the pier diameter, V is the flow velocity, V_c is the critical velocity, and Fr is the Froude number.

Table A. Well-known equations for the estimation of scour depth around bridge piers

No.	Reference	Equation
1	Inglis-Poona I (Inglis, 1949)	$d_{se}=1.7D^{0.22}{V}^{0.52}{Y}^{0.52}-Y$ All variables are in English units.
2	Inglis-Poona II (Inglis, 1949)	$d_{se}=1.73D^{0.22}{Y}^{0.78}-Y$
3	Laursen and Toch (1956)	${\displaystyle \frac{d_{se}}{Y}}=1.35{\left({\displaystyle \frac{D}{Y}}\right)}^{0.7}$
4	Chitale (1962)	${\displaystyle \frac{d_{se}}{Y}}=6.65Fr-0.51-5.49F{r}^2$
5	Blench-Inglis I (Blench et al., 1962)	$d_{se}=1.8D^{0.25}{Y}^{0.75}-Y$ All variables are in English units.
6	Blench-Inglis II (Blench et al., 1962)	$d_{se}=1.53D^{0.25}{V}^{0.5}{Y}^{0.5}{d}_{50}^{-0.125}-Y$ All variables are in English units except d50 which is in mm
7	Larras (1963)	$d_{se}=1.42K_s{D}^{0.75}$ Ks is pier shape factor All variables are in English units.
8	Neill (1964)	$d_{se}=1.5D^{0.7}{Y}^{0.3}$
9	Breusers (1965)	$d_{se}=1.4D$
10	Shen et al. (1969)	${\displaystyle \frac{d_{se}}{D}}=3.4(Fr{)}^{2/3}{\left({\displaystyle \frac{Y}{D}}\right)}^{1/3}$
11	Blench (1969)	${\displaystyle \frac{d_{se}+Y}{Y_r}}=1.8{\left({\displaystyle \frac{D}{Y_r}}\right)}^{0.25}$ $Y_r=1.48{\left({\displaystyle \frac{q^2}{F_B}}\right)}^{\frac{D}{3}}$ $F_B=1.9(d_{50}{)}^{0.5},d_{50}\text{in mm},qin\left({\displaystyle \frac{m^2}{s}}\right)$
12	Shen (1971)	$d_{se}=0.00022{\left({\displaystyle \frac{VD}{\nu }}\right)}^{0.619}$ $\nu$ is the kinematic viscosity of water.
13	Hancu (1971)	${\displaystyle \frac{d_{se}}{D}}=2.42\left(2{\displaystyle \frac{V}{V_c}}-1\right){\left({\displaystyle \frac{V_c^2}{gD}}\right)}^{\frac{1}{3}}$
14	Coleman (1971)	${\displaystyle \frac{V}{\sqrt{2g{d}_{se}}}}=0.6{\left({\displaystyle \frac{V}{D}}\right)}^{0.9}$
15	Neill (1973)	$d_{se}=K_{s}D$ Ks is pier shape factor.
16	Breusers et al. (1977)	${\displaystyle \frac{d_{se}}{D}}=f\left({\displaystyle \frac{V}{V_c}}\right)\left[2\tanh \left({\displaystyle \frac{Y}{D}}\right)\right]K_1{K}_2$ $f\left({\displaystyle \frac{V}{V_c}}\right)=\{\begin{array}{c}0{\displaystyle \frac{V}{V_c}}\le 0.5\\ \\ \left(2{\displaystyle \frac{V}{V_c}}-1\right)0.5<{\displaystyle \frac{V}{V_c}}<1\\ \\ 1{\displaystyle \frac{V}{V_c}}\ge 1\end{array}$ K1: pier shape factor K2: pier alignment factor
17	Jain and Fischer (1979)	${\displaystyle \frac{d_{se}}{D}}\{\begin{array}{c}2{\left(Fr-F{r}_c\right)}^{0.25}{\left({\displaystyle \frac{Y}{D}}\right)}^{0.5}\left(Fr-F{r}_c\right)\ge 0.2\\ 1.84F{r}_c^{0.25}{\left({\displaystyle \frac{Y}{D}}\right)}^{0.3}\left(Fr-F{r}_c\right)\le 0\\ max\left(2{\left(Fr-F{r}_c\right)}^{0.25}{\left({\displaystyle \frac{Y}{D}}\right)}^{0.5},1.84F{r}_c^{0.25}{\left({\displaystyle \frac{Y}{D}}\right)}^{0.3}\right)0<\left(Fr-F{r}_c\right)<0.2\end{array}$Frc is critical Froude number.
18	Froehlich (1988)	${\displaystyle \frac{d_{se}}{D}}=0.32K_{s}F{r}^{0.2}{\left({\displaystyle \frac{D_p}{D}}\right)}^{0.62}{\left({\displaystyle \frac{Y}{D}}\right)}^{0.46}{\left({\displaystyle \frac{D}{d_{50}}}\right)}^{0.08}$ Ks is pier shape factor, Dp is the projected width of pier.
19	Froehlich design equation (1988)	${\displaystyle \frac{d_{se}}{D}}=0.32K_{s}F{r}^{0.2}{\left({\displaystyle \frac{D_p}{D}}\right)}^{0.62}{\left({\displaystyle \frac{Y}{D}}\right)}^{0.46}{\left({\displaystyle \frac{D}{d_{50}}}\right)}^{0.08}+1$ Ks is pier shape factor, Dp is the projected width of pier.
20	Melville and Sutherland (1988)	${\displaystyle \frac{d_{se}}{D}}=K_d{K}_I{K}_y{K}_{\alpha }{K}_s$ Kd = sediment size coefficient; KI = flow intensity coefficient; Ky = flow shallowness coefficient; Ks = pier shape coefficient; Kα= pier alignment coefficient
21	Ansari and Qadar (1994)	$d_{se}=\{\begin{array}{c}0.024D^{3}D<7.2ft\\ 2.238D^{0.4}D>7.2ft\end{array}$
22	Simplified Chinese equation (Landers & Mueller, 1996)	$d_{se}=1.141K_s{D}^{0.6}{Y}^{0.15}{d}_{50}^{-0.07}\left({\displaystyle \frac{V-V_{ic}}{V_c-V_{ic}}}\right)$ $V_{ic}=0.645{\left({\displaystyle \frac{d_{50}}{D}}\right)}^{0.053}{V}_c$ All variables are in English units, and Ks is pier shape coefficient.
23	Melville (1997)	$d_{se}=K_{yD}{K}_I{K}_d{K}_s{K}_{\theta }$ KyD = depth-size coefficient; KI = flow intensity coefficient; Kd = sediment size coefficient; Ks = pier shape coefficient; K${}_{\theta }$ = pier alignment coefficient
24	Lyn (2008)	${\displaystyle \frac{d_{se}}{Y}}=\text{min}\left[0.21{\left({\displaystyle \frac{\text{V}}{{\text{V}}_{\text{c}}}}\right)}^{2.95},0.6\right]$
25	HEC-18 (Arneson et al., 2012)	${\displaystyle \frac{d_{se}}{Y}}=2K_1{K}_2{K}_3(Fr{)}^{0.43}{\left({\displaystyle \frac{D}{Y}}\right)}^{0.65}$ K1 = pier shape factor K2 = flow angle of attack factor K3 = bed condition factor
26	Sheppard et al. (2014)	$d_{se}=\{\begin{array}{c}2.5D{f}_1{f}_2{f}_{3}0.4\le {\displaystyle \frac{V}{V_c}}<1\\ D{f}_1\left[2.2\left({\displaystyle \frac{\frac{V}{V_c}-1}{\frac{V_{1p}}{V_c}-1}}\right)+2.5f_3\left({\displaystyle \frac{\frac{V_{1p}}{V_c}-{\displaystyle \frac{V}{V_c}}}{{\displaystyle \frac{V_{1p}}{V_c}}-1}}\right)\right]1\le {\displaystyle \frac{V}{V_c}}\le {\displaystyle \frac{V_{1p}}{V_c}}\\ 2.2D{f}_1{\displaystyle \frac{V}{V_c}}>{\displaystyle \frac{V_{1p}}{V_c}}\end{array}$ $V_{1p}=max\left(5V_c;0.6\sqrt{gY}\right)$ $f_1=f\left({\displaystyle \frac{Y}{D}}\right)$ $f_2=f\left({\displaystyle \frac{V}{V_c}}\right)$ $f_3=f\left({\displaystyle \frac{D}{d_{50}}}\right)$
27	Choi et al. (2017)	${\displaystyle \frac{d_{se}}{D}}=9.694K_s{\left({\displaystyle \frac{V}{\sqrt{gD}}}\right)}^{0.64}{\left({\displaystyle \frac{Y}{D}}\right)}^{0.056}{\left({\displaystyle \frac{d_{50}}{D}}\right)}^{0.082}{\left({\displaystyle \frac{V}{V_c}}\right)}^{1.53}$
28	Shamshirband et al. (2020)	${\displaystyle \frac{d_{se}}{Y}}=0.893(Fr{)}^{1.016}{\left({\displaystyle \frac{D}{Y}}\right)}^{0.625}{\left({\displaystyle \frac{d_{50}}{Y}}\right)}^{-0.262}{\left({\displaystyle \frac{V}{V_c}}\right)}^{-0.836}$