The dataset used in this study is the 2014 U.S. Geological Survey Pier-Scour Database (PSDb-2014), compiled by Benedict and Caldwell (2014) through a systematic literature review of published laboratory and field measurements of pier-scour. For the present analysis, the laboratory subset consisting of 569 data points was utilized. These data encompass a wide range of hydrologic, sediment, and geometric conditions, representing diverse flow regimes and experimental configurations. The dataset is widely regarded as a benchmark source for scour prediction studies and provides a reliable foundation for reproducible model development and comparison with prior research.

The equilibrium scour depth at a bridge pier is governed by three categories of parameters: (i) bed-material properties, (ii) water inflow conditions, and (iii) pier geometry (Pizarro et al., 2020). In alignment with this framework, the principal variables adopted in this study include the pier width normal to flow (b_n), approach flow velocity (V_o), sediment critical velocity (V_c), approach flow depth (y_o), median sediment size (D₅₀), geometric standard deviation of the sediment-size distribution (σ_g), and the duration of flow/scouring (T). Accordingly, the equilibrium scour depth at the pier (y_s) can be expressed as a function of these variables, as shown in Equation 1. y s = F b n , V o , V c , y o , D 50 , σ g , T (1)

The histograms in Figures 2A–H illustrate the distributional features summarized in Table 1, revealing the effective range of each variable. The pier width b n is tightly clustered at smaller sizes (0.05–3.00 ft; median 0.35 ft) with a strongly right-skewed tail. Approach and critical velocities ( V o and V c ) share similar shapes: most V o values lie below 4 ft/s but extend to 7 and 4.18 ft/s, with medians of 1.35 and 1.15 ft/s. Flow depth   y o Spans 0.07–6.23 ft (median 0.56 ft), concentrating in shallow flows. Median sediment size D 50 covers 0.22–7.80 mm (median 0.80 mm), heavily right-skewed toward fine–medium sands, while sediment gradation   σ g is narrowly distributed (1.10–1.50; median 1.30), indicating well-sorted beds. Duration T has the greatest range (29–54,780 min; median 2,998 min) with a strong tail. The target scour depth, y_s is right-skewed (0.01–4.63 ft; median 0.32 ft), with most events below 1 ft and a sparse upper tail. These distributional patterns in Figure 2 corroborate the summary statistics presented in Table 1.

Given the extensive positive skew across several variables ( b n , V o , V c , y o , D 50 , T , y s ) and heavy tails, particularly for T and y s , robust predictive models are recommended, as they tend to produce more stable and symmetric residuals (Yang et al., 2019). Although   σ g shows limited marginal variability, potential interactions with D 50 and flow variables may influence scour response. Because y o , V o , and D 50 span wide ranges (Table 1), models should be validated across hydraulic intensities and sediment coarseness (Baranwal and Das, 2024a). Dimensionless groups, flow intensity ( V o / V c ) and relative depth ( y o / b n ) are recommended for sensitivity analysis, as they can reduce scale effects and linearize correlations with y _s (Hassan and Jalal, 2021). Furthermore, outlier-aware evaluation metrics (e.g., median-based and quantile-based statistics) and uncertainty reporting are applied to prevent rare, long-duration tests or large-scour cases from biasing risk-relevant predictions (Arachchige and Prendergast, 2024).

The Pearson correlation heatmap (Figure 3) indicates two dominant patterns among the study variables. The critical velocity V c shows a strong positive correlation with median sediment size D 50 ( 0.95 ), and approach flow velocity V o relates to both ( 0.44 with V c and 0.43 with D 50 ), the same correlation was also reported by various studies (Chen F. et al., 2025; de Lange et al., 2024). By contrast, the pier width normal to flow bn has a strong correlation with y s ( 0.86 ) and a substantial relationship with time T ( 0.61 ), while T itself relates to y s   ( 0.66 ), validating the correlations reported by recent studies (Nandi and Das, 2025a; Xu et al., 2025). Initial flow depth y 0 shows a moderate relation with y s ( 0.53 ) and weaker ties to the hydraulics (V_o and V_c) block, having a correlation of 0.11 and 0.39. Sediment grain-size distribution ( σ g ) is largely independent of the other factors, slightly negative with V c and D 50 ( − 0.13 ) and negligible with y s ( − 0.12 ) (Baranwal and Das, 2024b). Overall, the correlation structure suggests that geometric configuration (bn and yo) and flow duration (T) are the primary predictors of scour depth, whereas velocity–size variables (V_o, V_c, and D₅₀) are strongly interrelated and should be treated carefully to mitigate multicollinearity in model development.

All variables were standardized with consistent symbols and units, including b_n (ft), V_o (ft/s), V_c (ft/s), y_o (ft), D₅₀ (mm), σg (−), T (min), and ys (ft) by trimming headers and mapping aliases; columns were then forced to numeric with placeholders treated as NaN (Jacobson et al., 2024; Kang, 2013; Peng et al., 2023). For the laboratory dataset, rows containing missing values were removed to ensure that descriptive statistics and model training were based exclusively on observed measurements. The resulting cleaned dataset was then stored and used to develop the scour-depth prediction models (Kang, 2013; Peng et al., 2023).

Because the compiled measurements originated from multiple laboratories, the raw distributions contain outliers and scale heterogeneity. We therefore applied normalization and standardization (feature scaling) before modeling to improve numerical conditioning and robustness. As evident in Table 1 and the histograms in Figures 2A–H, the dataset is biased; certain value ranges are over-represented while others are limited or absent, which can bias learning and increase the risk of overfitting (Charilaou and Battat, 2022). The adopted preprocessing workflow, including consistent labeling, type coercion, outlier-aware scaling, and missing-data handling, was therefore implemented to reduce these effects and enhance generalization performance in the predictive modeling stage.

For the Knik River piers, all records were retained, and units were harmonized across variables. A new hydraulic descriptor, the Froude number, was engineered as given in Equation 2 (CHOW, 1959). Fr = Vo / g  yo (2) where g = 32.174 ft s 2   . Remaining gaps in {b_n, V_o, y_o, D₅₀, σ_g, y_s, Fr} were imputed using distance-weighted k-nearest neighbors (KNN) imputation to preserve local structure within the feature space.

Two variables frequently missing from the records, V_c and T, were not statistically imputed. Instead, they were estimated using physics-integrated formulas within three defined hydraulic scenarios (baseline, worst-case, and extreme-quantile) to enable risk-aware model evaluation. For modeling, tree-based algorithms operated on raw-scale variables, whereas kernel-based algorithms applied in-pipeline standardization. All cleaned and scenario-augmented datasets were archived to ensure full transparency and reproducibility.

After preprocessing, the standardized laboratory dataset was randomly partitioned into 80% training and 20% testing subsets, following recommendations by Bichri et al. (2024). A fixed seed ensured full reproducibility. To better preserve the empirical distribution of the target variable (y_s), stratified random splitting was performed using binned scour depths, ensuring that both training and test sets captured the complete range of scour magnitudes and covariate combinations. This strategy mitigates selection bias and supports more reliable model generalization (Kapoor and Narayanan, 2023).

Model training and any preprocessing operations that could introduce data leakage (e.g., scaling for kernel-based methods) were encapsulated within scikit-learn pipelines and fitted exclusively on the training data. The held-out test data were never used during training, hyperparameter tuning, or preprocessing (Wu et al., 2025). This setup provides an unbiased assessment of out-of-sample performance while ensuring that both splits represent diverse hydraulic, geometric, and sediment conditions (Hameed et al., 2025).

Model development and analysis were performed in Google Colab (Python 3.12.12, Linux-6.6.105+), using NumPy 2.0.2, pandas 2.2.2, scikit-learn 1.6.1, SHAP 0.50.0, and Matplotlib 3.10.0. A total of five machine-learning models were developed and benchmarked, with their hyperparameters optimized using Random Search (RandomizedSearchCV, 50 random trials). The models were Gradient Boosting, AdaBoost (tree-based), XGBoost, Gaussian Process Regression with an RBF kernel, and Kernel Ridge Regression with a polynomial kernel. Each model was implemented within a leak-safe pipeline: tree ensembles (Gradient Boosting, AdaBoost with shallow trees, and XGBoost) used raw feature scales, whereas kernel-based methods (GPR with RBF and KRR with Poly) incorporated standardized features via StandardScaler to inputs. The five models, their algorithmic families, and principal hyperparameters are summarized in Table 2.

Model performance was assessed using three commonly adopted statistical indicators: the coefficient of determination (R²), root mean square error (RMSE) and mean absolute error (MAE) evaluated for both training and testing subsets to verify generalization and detect potential overfitting. These performance metrics are widely used in the literature for validating predictive models (Khajavi et al., 2025; Madurwar et al., 2025; Ramujee and Praseeda, 2025).

The coefficient of determination (R²) quantifies the proportion of variance in the observed data explained by the model, ranging from − ∞ to 1, where a value of 1 indicates a perfect fit and 0 corresponds to a mean-only baseline (Bentegri et al., 2025). RMSE and MAE are non-negative error measures in the same units as the target variable, with 0 representing perfect prediction (Mohammed et al., 2025). RMSE is generally equal to or greater than MAE, as it penalizes larger deviations more heavily.

Together, these three indicators provide a balanced evaluation framework: R² reflects the model’s explanatory power, RMSE captures overall prediction accuracy with sensitivity to large errors, and MAE represents typical error magnitude and robustness to outliers. The mathematical formulations for these evaluation metrics are presented below in Equations 3–5. R 2 = 1 − ∑ i = 1 n y i − y ^ i 2 ∑ i = 1 n y i − y – 2 (3) RMSE = 1 n ∑ i = 1 n y i − y ^ i 2 (4) MAE = 1 n ∑ i = 1 n y i − y ^ i (5) where y i is the observed value for sample i , y ^ i is the predicted value for sample i , y – is the mean of the observed values, n is the number of samples

K-Fold cross-validation partitions the available training data into K equally sized, non-overlapping folds and performs K repeated training and validation cycles. In each cycle, one-fold serves as the validation set while the remaining (K-1) folds constitute the temporary training subset (Teodorescu and Obreja Braşoveanu, 2025; Wilimitis and Walsh, 2023). This rotation ensures that every observation is used for validation exactly once and for training (K-1) times, producing a distribution of performance scores that reflects model sensitivity to data variability and the bias–variance trade-off (Kapoor and Narayanan, 2023).

To estimate out-of-sample performance and minimize split-specific bias, 5-fold shuffled cross-validation was applied on the training set only (representing 80% of the total data), consistent with the approach used by Al-Shamasneh et al. (2025). The training data were stratified by scour depth ranges so that each fold maintained a similar distribution of the response variable, ensuring balanced representation across folds. A fixed random seed was employed to ensure repeatability.

For each fold, models were trained on 80% of the training data and validated on the remaining 20%. All preprocessing steps, including scaling for kernel-based methods, were encapsulated within the scikit-learn pipelines to prevent data leakage. Performance was recorded as R² and RMSE for each fold, allowing for a consistent comparison of model stability and predictive accuracy across the five folds.

Sensitivity analysis examines how variations in model inputs influence the predicted output. In this study, it was applied to quantify how changes in key hydraulic, geometric, and sediment variables affect the predicted equilibrium scour depth (y _s ). The analysis employed SHAP (SHapley Additive ExPlanations), a principled framework for both global and local model interpretability based on cooperative game theory.

SHAP treats each input feature as a “player” in a game whose contribution to the prediction is computed as the average of its marginal effects across all possible future combinations. The resulting Shapley values satisfy the desirable properties of local accuracy, missingness, and consistency, making SHAP a unique and theoretically grounded additive explanation model (Aas et al., 2021). This framework unifies numerous prior feature-attribution methods and provides both local explanations (case-specific contributions that increase or decrease predictions) and global explanations through mean absolute SHAP values aggregated across samples.

TreeSHAP extends this formulation to decision trees and ensembles, enabling polynomial-time computation of exact Shapley values while also accounting for feature interaction effects. It provides practical visualization tools such as mean SHAP bar charts, beeswarm summaries, and dependence plots that capture local-to-global behavior without compromising fidelity to model predictions (Lundberg et al., 2020; Lundberg and Lee, 2018).

Overall, SHAP provides a theoretically sound and practically effective bridge between per-instance reasoning and global sensitivity of complex models. When applied with consideration for inter-variable dependence and supported by complementary diagnostics, it enables transparent, reproducible insights into how hydraulic, geometric, and sediment parameters collectively influence predicted scour (Alasmari et al., 2025).

Bridge scour arises from the interaction of hydraulic intensity, event duration, sediment mobility, and pier geometry. Field datasets rarely span the full range of conditions that govern safety, and key drivers such as the critical velocity for bed mobility (V_c) and the effective duration of mobility (T) are often missing (Belmokhtar et al., 2025; Shanmugam et al., 2025). A scenario framework allows us to impute defensible values where observations are incomplete by implementing physics-informed integration into the models. The framework further propagates credible extremes through a physics-guided predictor of scour depth (y_s) and summarizes risk per pier by taking the envelope across scenarios. Scenarios convert patchy records into decision-ready evidence about plausible and upper-bound scour. To operationalize the analysis, this study considers four application scenarios, Q99, WC-VcT, WC-Flow + Base, and WC-Flow + VcT.

The Knik River bridge in southcentral Alaska is examined as a case study to estimate credible upper bounds on pier scour and to cross-validate structural sufficiency under forecast and design-flood conditions. The proposed methodology estimates extreme pier scour under rare and persistent flood conditions by developing four physically consistent, pier-specific scenarios: Q99, WC-Flow, WC-VcT, and WC-Flow + VcT. Scour depth is predicted after recalculating the mobility (sediment motion) threshold velocity and event duration to reflect changes in hydraulic and sediment states. Using this setup, we can fold V_c and T into the scenario envelope via physics-informed learning model integration and validate them in cases with missing V_c or T values. This technique shows the model can both predict pier scour depth and impute the missing parameters in a physics-consistent way. The process begins with standardized inputs for each pier, including b n , V o , V c , y_o, D₅₀, σ g , T, and y_s. All data are converted to numeric values, with placeholders treated as missing and filled using physics-based formulas before the scenario modifications are applied.

Values of V c is estimated using a Shields-type formula that considers sediment gradation (Soulsby, 1997; van Rijn, 1984). The formula for calculating the critical velocity can be expressed as Equation 6. V c = C uc g · D 50 , ft · σ g α (6)

Where C uc = 2.6 , α = 0.15 , g = 32.174   ft / s 2 ,

and D 50 , ft = D 50 , mm × 10 − 3 × 3.28084

Values of duration of flow (T) are estimated from an advection-based time scale modified by the mobility ratio, using T = 60 y o V o V o V c β with β = 1.25 , ensuring that duration stays between 20 min and 3 weeks (Melville and Chiew, 1999). For other inputs ( b n , V o , y o , D 50 , σ g ), gaps are filled using robust median values. To avoid uniformity across the dataset, each pier’s rank within the empirical distributions of V o , y o , D 50 , and σ g is calculated so that high-energy piers receive stronger scenario variation. After any variable changes, both V c and T are recalculated to ensure consistency with the modified state. The updated threshold velocity is expressed in Equation 7 (Julien, 2010; Soulsby, 1997; van Rijn, 1984). V c ← V c , nominal y o y o 0 − 0.12 (7) where V c , nominal = C uc g · D 50 , ft σ g α .

The event duration is adjusted accordingly using T ← 60 y o V o V o V c 1.35 y o y o 0 0.55 . These relationships make the recalculated values dynamically consistent with the altered hydraulic and sediment conditions. The Q99 scenario increases each pier’s V o , y o , D 50 , and σ g toward their 99th-percentile values, scaled by pier rank so that higher-energy piers experience stronger effects. This represents rare but plausible extremes while maintaining pier-level diversity. In the WC-Flow scenario, both V o and y o are multiplied by rank-dependent factors with small random noise, as given in Equations (8) and (9) (Brandimarte et al., 2012; Das et al., 2021). V o ← V o 1.45 + 0.60 r V o 1 + 0.03 ε  and (8) y o ← y o 1.35 + 0.45 r y o 1 + 0.02 ε (9)

Sediment characteristics are forced toward their upper distribution tails, and V c and T are then recalculated to capture strong, short-duration hydrodynamic forcing. The WC-VcT scenario focuses on situations in which sediment motion is sustained through long-lasting floods or multiple storm events. It lowers V c and increases T based on pier rank, as given in Equations 10,11 . V c ← V c 0.65 − 0.15 r V o  and (10) T ← T 3.0 + 3.0 r V o (11)

While applying smaller adjustments to V o and y o to maintain variability across piers. The combined WC-Flow + VcT scenario merges both sets of conditions, applying the strong hydrodynamic changes from WC-Flow and the persistence adjustments from WC-VcT. After applying these combined multipliers to V o , y o , V c , and T , both values are recalculated again to maintain physical consistency, producing the most conservative condition for resilience checks. Scour depth for each scenario is estimated using an ensemble of Gradient Boosting, AdaBoost (Tree), XGBoost, Gaussian Process (RBF), and Kernel Ridge (Polynomial) models. The input vector is x = b n , V o , V c , y o , D 50 , σ g , T . The physics-informed formulation used to predict pier scour depth is given by Equation 12 (Bruce and Melville, 2000). y ^ s = k   b n   Fr m   y o n (12) with Fr = V o g y o , m = 0.38 , and n = 0.33

For each pier, the maximum estimated scour depth across all scenarios, y s max = max ⁡ scenario y ^ s , is selected and paired with the scenario that produces it. In interpretation, the Q99 scenario provides hazard-based checks for extreme statistical conditions, while WC-Flow represents short, intense flood events. The WC-VcT scenario highlights the effects of prolonged sediment mobility under persistent floods, and WC-Flow + VcT defines the upper bound useful for countermeasure design and durability assessments. Since V c and T are recomputed for each scenario using Equations 13, 14, the model reflects condition-specific hydraulics. V c ∝ g   D 50   σ g α y o / y o 0 − 0.12 and (13) T ∝ y o V o V o V c 1.35 y o y o 0 0.55 (14)

By applying this methodology, all scenarios remain physically realistic. The resulting scour envelopes reveal whether the peak intensity of flow or prolonged duration dominates pier scour risk for each pier of the Knik River bridge. The resulting pier-specific values are summarized in Table 3 (A), V_c (ft/s) by pier across each scenario, and Table 3 (B), T (min) by pier across each Scenario. Together, these tables show how mobility thresholds drop, and durations lengthen as scenarios intensify. These patterns directly amplify predicted scour and help flag piers most vulnerable under extreme, long-lasting floods.

The developed models were evaluated based on the performance indicators, including R², RMSE, and MAE, which were also implemented by other researchers for ML model performance evaluation (Khoshvaght et al., 2025; Koçak, 2025; Shobayo et al., 2025). A higher R² (closer to 1) means the model explains more of the variation and fits the data better (Mamudu et al., 2025). Lower RMSE and MAE reflect higher accuracy, and values of 0 correspond to a perfectly fitting model (Chai and Draxler, 2014).

Across models, performance is uniformly strong, with tree ensembles leading as illustrated in Figures 4A–E plots. Gradient Boosting delivers the top generalization, having the best fit line R² (0.98), and the parity plot shows predictions clustered tightly around the 45° line, as shown in Figure 4A. The best-fit line nearly overlaps the parity line over most of the range, indicating minimal systematic bias, with only a few high-value points pulling slightly above the diagonal. AdaBoost (Tree) closely mirrors this behavior, having fit-line R² = 0.97, showing a similarly tidy spread and a best-fit line that tracks the diagonal with small deviations at the upper end, as depicted in Figure 4B. The developed XGBoost model parity plot is largely well aligned. The best-fit line for XGBoost remains close to the 1:1 line but is influenced by a handful of larger residuals at high targets, which explains the higher test RMSE. Kernel methods provide stable baselines with smooth behavior. Gaussian Process (RBF) and Kernel Ridge (Poly), have fit line R-squared value of 0.97 and 0.95, maintain compact clouds near the origin and a gradual, orderly dispersion as values grow. In both plots, the best-fit line sits just below the parity line at higher observed values, reflecting mild underprediction in the extreme range while remaining well calibrated through the bulk of the data.

Moreover, to assess overfitting and underfitting, the gap between training and testing metrics should be small across all performance indicators (Aliferis and Simon, 2024; Emmert-Streib and Dehmer, 2019). Table 4 summarizes the training and testing performance of the developed models. Gradient Boosting and AdaBoost (Tree) exhibit small training-testing gaps, high testing R² (0.96) with modest increases in RMSE and MAE, indicating good generalization. XGBoost shows slight signs of overfitting, with near-perfect training metrics having R² (0.999) and deteriorating more on the testing set R² (0.939), higher RMSE, and MAE. The resulting difference in R² between training and testing is 0.06, indicating that the model performs very well overall, with only a small number of test samples acting as outliers that fall noticeably below the 1:1 parity line. The kernel methods, Gaussian Process (RBF) and Kernel Ridge (Poly), yield the largest test errors and lower test R² (0.927–0.929), reflecting a slightly weaker generalization. Overall, the parity visuals and summary table agree that the predictions are well centered with tight residual structure for most targets, and the best-fit lines lie closely to the ideal diagonal across models, most notably for Gradient Boosting and AdaBoost, highlighting a strong capture of all the data points.

Residual plots (observed minus predicted) are used to screen for structural non-linearity, outliers, and to compare behavior on training versus testing data as a check on generalization (Kumar et al., 2025; Sharma et al., 2025). Checking residuals on the test and training sets helps separate true signal from overfitting, and models that generalize well show small, random-looking residuals in both. Figures 5A–E depicts residual histograms (observed vs. predicted) for all five models are tightly centered near zero, with substantial overlap between the training (blue) and testing (red) curves, indicating that none of the models is markedly overfitting. Gradient Boosting shows a compact, symmetric spread around zero, suggesting low bias and stable variance. AdaBoost (Tree) and Gaussian Process (RBF) exhibit a mild right-shift of the test curve (slightly positive bias), but the displacement is small relative to their overall dispersion. XGBoost produces a narrow core with light right tails, implying good central accuracy with a few underpredicted cases. Kernel Ridge (Poly) has slightly broader tails than the boosting models, though its training–testing overlap remains high. Overall, the distributions are unimodal and roughly symmetric, residual magnitudes are modest, and the similarity of train/test shapes supports good generalization across the range of scour depths represented in the data.

For the Gradient Boosting model, residuals r = y obs − y pred were used to derive empirical prediction intervals. The 5th and 95th percentiles of the residuals are − 0.0449 ft and 0.0401 ft on the training set and − 0.0949 ft and 0.1723 ft on the testing set, respectively. These values imply an approximate 90% prediction interval of y ^ s + − 0.095 , 0.172 ft on unseen data, corresponding to a typical uncertainty band of about ± 0.13 ft around each predicted scour depth. Together with the close agreement between train- and test-based bands, this narrow interval width indicates that the Gradient Boosting model achieves accurate and well-calibrated scour depth predictions with only modest uncertainty on unseen cases.

The actual vs. predicted plots, as shown in Figures 6A–E, demonstrate how well each model tracks measured scour depths across the whole sample index, also utilized by other studies (Nandi and Das, 2025b; Showkat et al., 2025). Overall, alignment is strong for Gradient Boosting and AdaBoost, whose test traces are close to the observed values with very minor variations in the top tail, indicating good generalization from low to high values. XGBoost follows the pattern well but demonstrates the previous overfitting tendency; training points are practically precise, but a few significant (y_s) test spikes are slightly under-predicted. The kernel approaches, Gaussian Process (RBF) and Kernel Ridge (Poly), recreate the central range stably while smoothing the highest peaks more than the boosting models, which is consistent with their slightly higher test errors. Errors are concentrated at the largest (y_s) events across all the models, as is predicted in hydraulics, where extremes are infrequent and difficult to learn, but most mid-range depths are well captured (McInerney et al., 2020).

The current study utilizes K-fold for cross-validation to estimate out-of-sample performance by splitting the data into k folds, training on k − 1 folds, and evaluating on the held-out fold, then repeating so each sample is tested once and averaging the scores. The K-Fold cross-validation is also employed by various researchers to validate their ML models performance against overfitting or underfitting, indicating its usefulness and importance in the ML models performance evaluation (Habal and Benbouras, 2025; Jamal and Ahmed, 2025; Zhao et al., 2025). As shown in Table 5, ensemble-based models perform better and more consistently when cross-validated. Gradient Boosting and XGBoost achieve higher (R²) while maintaining low, densely dispersed RMSE, with XGBoost having the greatest single-fold result (R² = 0.95) and Gradient Boosting demonstrating somewhat greater stability across five folds (R² = 0.83–0.94). AdaBoost is competitive but has more variability. The kernel models are sensitive to how the data are split. Gaussian Process varies a lot between folds and shows the clearest outliers, including one-fold with low (R²) and high RMSE. Kernel Ridge has some strong folds, but overall, it is less consistent and shows a few moderate outliers. On the other hand, the tree ensembles show compact error distributions; XGBoost and AdaBoost each exhibit an isolated high-error fold, whereas Gradient Boosting maintains the tightest spread with only minor deviations. Overall, Table 5 validates that ensemble tree approaches provide the most robust generalization on this dataset, while GP and KR are comparatively volatile.

To identify the global drivers of pier-scour depth, SHAP (Tree SHAP) was applied to the trained models to decompose predictions into additive feature contributions and aggregate them across the dataset. This yields a ranked importance profile, supported by beeswarm and dependence plots that highlight which hydraulic, geometric, and sediment variables most consistently increase or decrease predicted scour. SHAP is now widely used to quantify global feature importance and has become a leading tool for model interpretation, with recent studies showing it provides consistent, dataset-wide insights into how predictors drive outputs (Cappelli et al., 2023; Cappelli and Grimaldi, 2023; Mushtaq et al., 2024).

Figures 7A, B depicts SHAP feature-importance and SHAP beeswarm, both of which jointly demonstrate that the pier width normal to flow b_n (ft) is by far the dominant predictor of scour depth y_s (ft), accounting for 70.6% of the model’s overall importance. Flow speed V_o (ft/s) with 9.8% and approach depth y_o (ft) with 7.5% offer the next highest contributions, followed by event duration T (min) with 4.3%, and sediment gradation (σ_g) with 4.0% respectively. On the other hand, critical velocity (V_c) and sediment size (D₅₀) have a minor global influence of 2.3% and 1.4%. Moreover, the beeswarm plot, as shown in Figure 7B, clarifies directionality and nonlinearity at the observation level. Large bn values (warm points to the right) consistently increase predicted y_s, while small bn values reduce it, an effect that is strong and monotonic, also reported by Baranwal and Das (2024a), Fuladipanah et al. (2023). Also, higher V_o tends to shift predictions upward, while lower V_o tends to shift them downward, aligning with standards that faster approach flow promotes scour. The SHAP plot shows the effect of y_o is more moderate and slightly nonlinear; higher depths generally push ys upward, but with a visible spread that suggests interactions with V_o and b_n. The plot also depicts that longer T shows a mild positive trend (more exposure leads to more scour growth), also validated by Melville and Chiew (1999). Sediment-size distribution (σ_g) exhibits mixed local effects (both signs), consistent with gradation influencing scour depth, aligning well with Mir et al. (2018) conclusions. In contrast, higher V_c typically reduces predicted scour (points with high Vc cluster on the negative SHAP side), reflecting that a bed that requires larger velocities to mobilize is less prone to scour under the same forcing, as reported by Arneson et al. (2012). D₅₀ effects are small and mostly negative in this sample, indicating limited incremental predictive power once bn, V_o, and y_o are known. Taken together, Figures 7A, B) indicates a physically consistent hierarchy, geometry (b_n) dominates, hydraulics (V_o, y_o, T) provide substantial but secondary control, and mobility metrics (σ_g, V_c, D₅₀) modulate scour at the margins. The tight, mostly one-sided SHAP pattern for b_n and V_o also suggests the model learned stable, interpretable relationships rather than relying on incorrect interactions.

To evaluate the consistency between classical correlation analysis and interpretable feature importance, a feature-correlation versus SHAP comparison plot is given in Figure 7C. The horizontal axis reports the Pearson correlation between each input and the target (y_s), and the vertical axis reports the global SHAP importance as a percentage of the total contribution. Overall, the two measures are consistent, with (b_n) exhibiting both the highest correlation and the largest SHAP importance, confirming its dominant influence on the model predictions. Variables such as (y_o) and (T) show moderate positive correlations and intermediate SHAP contributions, whereas weakly correlated inputs such as (D₅₀) and (V_c) contribute only marginally. The gradation parameter (σ_g) displays a slightly negative correlation with (y_s) but only modest SHAP importance, indicating that its effect is small and predominantly inverse. Overall, this comparison indicates that the Gradient Boosting model’s learned importance structure is broadly aligned with the underlying statistical relationships in the data.

Figures 8A–D presents SHAP dependence plots, where the x-axis corresponds to the value of the input feature and the y-axis to its SHAP value, representing feature’s marginal contribution to the predicted scour depth (y_s). Point colors encode the value of a secondary feature (as indicated by the accompanying color bar), with cooler tones represent lower values and warmer tones higher values. Systematic color gradients that coincide with changes in SHAP values highlight potential interaction effects between the two features. Within these plots, pier width (b_n) is the dominant predictor. SHAP values increase almost monotonically with (b_n) as shown in Figure 8A, implying larger expected scour for wider piers (Fuladipanah et al., 2023). Approach flow velocity (V_o) exerts a strong positive effect that saturates at higher speeds. Whereas initial depth (y_o) contributes positively but with greater dispersion, consistent with a secondary role modulated by geometry and hydraulics, as evident from Dong et al. (2025) work. Event duration (T) shows a threshold response, rapid SHAP increases from short to moderate durations, followed by a plateau, indicating diminishing marginal effects for long exposures. The colour overlay clarifies interactions: higher (V_c) and larger (b_n) elevate the SHAP contributions of (V_o) and (y_o), while smaller (D₅₀) (finer sediment) aligns with larger SHAP values for (T), meaning duration is more consequential on easily mobilized beds. Together, these patterns support the importance order (b_n > V_o > y_o > T) and reveal nonlinear, partially saturating responses shaped jointly by geometry, hydraulics, and sediment properties.

This section evaluates the generated scenarios, comparing predicted scour responses across stress-tested combinations of hydraulic, geometric, and sediment conditions to identify sensitivity patterns, potential worst cases, and the robustness of model conclusions. The predicted scour depth at Knik River bridge piers across all scenarios is presented in Table 6, with the frequency histograms in Figure 9A showing a clear ordering across scenarios. Among all the scenarios, the combined worst case (WC-Flow + VcT) dominates the upper range with a right-shifted distribution and a long upper tail, indicating both higher typical scour and more frequent extremes (several peaks around 4.7–4.9 ft). The WC-Flow scenario is still severe but lies slightly below the combined case, with a distribution centered at smaller depths and a shorter upper tail, consistent with peak hydrodynamic forcing without the extra scour development from extended duration. The WC-VcT scenario concentrates at lower depths with less spread, meaning longer periods above the mobility threshold increase scour, but the increases are smaller than those produced by peak-flow events. The Q99 scenario clusters narrowly around 4.0–4.3 ft and typically lies 0.3–0.6 ft below WC-Flow + VcT, marking a rare-but-plausible benchmark distinct from the engineered worst case.

The summary statistics plot, as shown in Figure 9B illustrates the distributions of scour-depth across the scenarios. Scenario WC-Flow + VcT yields the highest medians (3.8–4.3 ft) and the widest interquartile (IQR) ranges (1.2–1.6 ft), confirming both elevated typical scour and variability. Scenario WC-Flow shows slightly lower medians (3.4–3.9 ft) and narrower IQRs (0.9–1.3 ft). In many cases, medians differ from the combined scenario by only 0.1–0.3 ft, indicating sites where peak velocity is the primary driver. On the other hand, scenario WC-VcT produces distinctly lower medians (2.0–2.5 ft) with tighter IQRs (0.6–0.9 ft), highlighting a milder central tendency despite the role of duration of event (T). Scenario Q99 remains high but stable, reinforcing its use as a decision threshold separating climatological extremes from design-envelope stress tests.

The per-pier comparison depicted in Figure 9C reveals heterogeneous sensitivity. At several piers, WC-Flow nearly matches WC-Flow + VcT (differences lie 0.1–0.3 ft), signaling locations where peak flow alone governs risk and thus may warrant rapid-response triggers tied to rising velocity. At some piers, the larger gap between WC-Flow + VcT and WC-Flow signals a strong duration (V_c, T) effect. In other words, these sites are therefore more vulnerable to extended floods or multi-storm sequences and warrant long-term monitoring and additional protective measures. In contrast, at sites where WC-VcT matches or exceeds the other scenario, it indicates that what matters most is how long the strong flow lasts, not just how high the single peak flow. In other words, longer periods of strong current can dig more scour even if the peak is not the biggest, so watch how long the current stays strong enough to move sediment. Therefore, tracking how long the flow remains above the sediment-mobility threshold is essential, and this can be operationalized using the scenario envelope.

Across the figures, WC-Flow + VcT consistently produces the largest and most variable scour depths, WC-Flow is a close second at many piers, WC-VcT is lower but still consequential for duration-sensitive sites, and Q99 offers a compact, high-end benchmark below the engineered worst case. Together, these results support a tiered risk triage: scale design envelopes to the combined worst case, deploy rapid-trigger monitoring where peak flow dominates, and prioritize duration-aware mitigation where persistence drives risk. These results show that the proposed framework remains accurate and interpretable across diverse hydraulic conditions and scenario stress tests.

To facilitate practical implementation, an interactive bridge pier scour prediction application was developed and made available on the Hugging Face platform. The tool combines five calibrated models: Gradient Boosting, AdaBoost, XGBoost, Gaussian Process Regression, and Kernel Ridge Regression, in a unified interface. Users enter seven required input variables, including pier width normal to flow (b_n), approach flow velocity (V_o), sediment critical velocity (V_c), approach flow depth (y_o), median sediment size (D₅₀), geometric standard deviation of the sediment-size distribution (σ_g), and the duration of flow/scouring (T), all constrained to their allowable ranges. For each set of input values, the application provides the predicted scour depth from every model, allowing direct comparison of model outputs. All trained model files as well as the Python scripts for model training, evaluation, and SHAP-based interpretation are included in the Hugging Face repository. This structure promotes accessibility for practitioners, guarantees transparency and reproducibility for researchers, and is fully usable without any prior knowledge of coding. The link to the developed application, model scripts, and SHAP analysis code is provided in the data availability section of the manuscript.

Table 7 provides a quantitative comparison of the present models with established bridge-pier scour predictors spanning various methodological scopes. Among the developed models, the tree-based algorithms, the Gradient Boosting model achieves superior generalization, evidenced by the highest testing R² value (0.959) and the lowest test RMSE (0.145). AdaBoost offers closely comparable performance (test R² 0.958, RMSE 0.146). Relative to cylindrical-pier models reported in Fuladipanah et al. (2023), including MARS, GEP, and M5 model tree, the present tree-based regressors demonstrate clear improvement in both test accuracy and error rates (testing R² for MARS, GEP, and M5 were 0.917, 0.872, and 0.698; RMSE were 0.090, 0.114, and 0.284, respectively). Similarly, complex-pier models such as those developed by Tien Bui et al. (2020) achieved the highest accuracy of R²= 0.91 in training for ANN among the developed models and dropped to 0.82 in testing. Overall, the current boosted tree approaches yield a more balanced combination of predictive performance and model simplicity for predicting bridge-pier scour.

The SHAP analysis in this study reveals a clear hierarchy among predictors. Geometry, especially the width parameter b n , has the largest, consistent influence on scour depth. Hydraulic factors like approach flow velocity ( V o ) and depth ( y o ) also have positive impacts but show some nonlinearity. Mobility and exposure variables, such as duration ( T ), sediment gradation ( σ g ), critical velocity ( V c ), and sediment particle size ( D 50 ), have smaller or mixed effects. These results align with some recent studies: Eini et al. (2023) found pier diameter ( D ) most important, with flow depth ( Y ) and velocity ( V ) secondary, and sediment size ( d 50 ) and critical velocity negatively correlated. Emami et al. (2025) highlighted dimensionless time ( Ut / D ) and relative pier size ( D / D m ) as key for time-dependent scour. Large-scale analyses also identify shape factor ( K 1 ) and depth ratio ( y / b ) as top predictors, confirming the geometry-led control with hydraulics second and mobility minor (Piraei et al., 2025). While interpretations broadly align with prior work, the present study explicitly quantifies the effects of seven input features, whereas some studies restrict analysis to a minimal set, risking the exclusion of an influential variable whose true effect may be larger.

A unique feature of the present study is the application of a comprehensive scenario envelope (Q99, WC-Flow, WC-VcT, WC-Flow + VcT) for stress testing, which supplies robust, decision-ready outputs that are not present in previous related work. Taken together, the ML techniques high out-of-sample accuracy, minimal training–testing performance gaps, and scenario-driven stress testing set this approach apart from prior studies, while aligning outcomes with core hydraulic principles and practical engineering expectations. Moreover, the present study developed an interactive practical tool for bridge pier–scour prediction, enabling direct use of the trained models by practitioners. As summarized in Table 7, the existing bridge-scour literature does not offer a comparable implementation-oriented tool alongside its modeling frameworks.

Scour depth is the lowering of the riverbed around a bridge pier and is a key parameter for bridge resilience. During severe floods, the stage and velocity rise quickly, forming intense downflow, horseshoe, and wake vortices around the pier. These vortices concentrate shear, mobilize the surrounding sediments, and carry them away, exposing foundations. If the scour depth grows beyond design limits, the pier’s capacity is reduced, and structural failure can occur (Lee and Hong, 2019). This makes scour depth central to flood-resilient design and operations. It must be monitored in time, with procedures to track scour growth and issue early warnings when thresholds are approached. Historically, scour has been the leading cause of bridge failures, with flood-driven high stages and velocities acting as the primary driver of scour development, explored by the current study. Floods cannot be prevented, but their impacts can be managed. At the design stage, prediction of scour depth for the site using hydraulic conditions and credible flood forecasts is crucial. While in operation, the identification of piers at higher risk, scheduling timely strengthening or protection, and, when necessary, temporarily removing a bridge from service to avoid catastrophic outcomes is the need of the day.

The current study goes beyond building a single predictor of pier-scour depth and undertakes a full model development, validation, and interpretation cycle aimed at practical deployment. Multiple machine learners were trained and tuned, and their performance was quantified on both training and testing splits using complementary metrics (R², RMSE, MAE) to expose accuracy, bias, and dispersion (Rana et al., 2025). To guard against optimistic estimates tied to a particular split, the study ran K-fold cross-validation, summarizing fold-wise scores and variability. The cross-validation framework revealed the metrics to diagnose overfitting (high train/low test, large fold variance) and underfitting (uniformly low scores) rather than relying on a single, whole-dataset score (White and Power, 2023). Residual and parity plots were reviewed alongside the metrics to verify that errors were pattern-free, and that upper-tail behavior was understood. Also, the study emphasized model insights, not only predictions. The employed framework aggregated explanations across the fitted models and generated SHAP analyses to rank the global importance of the hydraulic, geometric, and sediment variables for predicting scour depth and probe local dependence and interactions, e.g., how the effect of approach velocity changes with pier width or initial depth for a specific pier or duration (Nandi and Das, 2025a). These explanations link the learned relationships to hydraulics, help identify regime-dependent behavior (peak-dominated vs. duration-sensitive), and provide actionable levers for design and operations.

Moreover, the workflow was designed for field use. Starting from curated USGS records, targeted feature engineering is applied (including physics-informed estimates of missing drivers such as V_c and T, trained and validated models with transparent checks, and then wrapped the results in a scenario envelope that supports forecasting, triage of at-risk piers, and targeted monitoring or reinforcement. The emphasis is on reproducible, decision-ready outputs rather than record-keeping alone. Furthermore, the study outlines a forward path: incorporating additional toolkits (uncertainty quantification, conformal prediction, dynamic (V_o/V_c) exposure metrics), expanding site-specific scenarios as data improve, and deepening physics–ML integration to better represent near-threshold mobility and extreme events. Taken together, the approach demonstrates how machine learning and AI can be applied systematically and responsibly to strengthen flood-resilient bridge design and operations, while leaving clear hooks for future refinement.

To check whether the hydraulic structures remain adequate under forecast and design flood conditions, the current study used Knik River bridge piers located in southcentral Alaska as a case study to establish conservative upper limits on pier scour. Scenario analysis translates predictive insights into actionable strategies for flood-resilient bridge design (Kosič et al., 2023). Each pier is tested against four critical stress scenarios: Q99 (representing rare high extremes of V_o, y_o, D₅₀, σ_g), WC-Flow (short, intense floods), WC-VcT (events with low V_c and long T), and WC-Flow + VcT (the most conservative envelope). For each scenario, input variations trigger re-computation of V_c and T, ensuring dynamic consistency in predictions. Design recommendations typically emerge in two forms. Sites where peak flows dominate are best managed with measures that dissipate velocity or split flow, such as guide banks and pier upgrades. While sites sensitive to duration require countermeasures focused on resisting prolonged flood exposure, like strengthening the piers, deeper cutoffs, and toe protection (Lagasse et al., 2001). The WC-Flow + VcT scenario provides an upper benchmark for intervention: piers with predicted y_s above established safety limits under this scenario should be prioritized for structural upgrades or heightened flood monitoring. This risk-based framework allows authorities and practitioners to strategically allocate engineering resources by identifying the governing scour scenario for each pier. Sites vulnerable to peak flows and sites sensitive to event duration can be managed with scenario-tailored toolkits, optimizing expenditures and maximizing risk reduction. Importantly, over-design can be avoided where existing resilience is sufficient, as the scenario analysis clarifies which piers require immediate strengthening and which can be safely monitored without strengthening. This targeted approach provides a robust pathway for enhancing bridge safety and longevity, especially under future flood uncertainty.

The study also uses a physics-informed scenario envelope to handle missing inputs and to turn forecasts into actionable risk. When key variables are unavailable at a site, most often the critical velocity (V_c) and the effective duration of mobility (T) can be computed by establishing hydraulic relations and embedding those formulas inside the scenario framework. By doing so, each scenario (e.g., peak-flow, duration-focused, combined worst case, or Q99) carries internally consistent approach flow velocity (V_o), (y_o), (V_c), and (T), allowing the model to estimate scour even where records are incomplete. Using forecast flood stage and velocity, the scenario envelope identifies which piers will face high intensity (V_o/V_c) and long exposure (large T), ranks them by predicted scour, and triggers early warnings for high-risk assets.

The approach is transparent, safe, and sustainable because it ties imputation to physics rather than ad hoc guesswork, remains usable when monitoring gaps exist, and scales naturally to real-time operations. In this way, the study contributes to flood resilience by combining publicly available hydrologic data, physics-based estimation of missing drivers, and scenario-based analytics to guide timely protection, monitoring, and communication before hazardous conditions develop.

The models are based solely on laboratory data, and while laboratory experiments cover many scenarios, caution is needed when applying these results to field conditions. Factors like scale effects, armoring, debris or ice impacts, and live-bed sediment movement are less represented in lab settings compared to real-world environments. For the Knik River site, certain input variables, especially V_c and T, were either estimated or set according to scenario definitions, which introduces uncertainty to the predictions. The Q99 scenario draws on percentile ranks within the dataset, not on a hydrological basis; more realistic predictions could be achieved by linking these scenarios to basin-specific flood frequency and sediment mobility models. Lastly, the relatively limited variation in D₅₀ and σ_g data restricts the model’s ability to fully explore the influence of sediment gradation.

Field calibration and validation at bridges with high-quality monitoring should be prioritized, along with Bayesian or bootstrap uncertainty quantification for y s envelopes. Physics-guided ML (e.g., constraints from threshold mobility and time-to-equilibrium) can reduce extrapolation risk in sparsely sampled regimes. Replacing ad hoc extremes with probabilistic scenario generation tied to joint flood-sediment statistics would align envelopes with target reliability levels for flood-resilient design. Finally, integrating per-pier envelopes and uncertainty into risk-cost decision tools would help agencies prioritize mitigation and monitoring where it yields the greatest resilience gains.