## 2 Materials and Methods  
An integrated methodology is employed in this study, comprising two complementary components. First, advanced three-dimensional CFD simulations are used to model the complex fluid-sediment interactions under a range of operational parameters. Second, a robust multivariate polynomial regression analysis is conducted to systematically quantify the influence of these parameters on sediment resuspension. FLOW-3D Hydro software was used in this analysis to build a 3D hydro-morphodynamic numerical model. Flow Science, Inc. created FLOW-3D Hydro, a commercial software kit. Flow-3D uses a finite volume approach to solve the continuity equation (Equation (1)) and the unsteady Reynolds-averaged Navier-Stokes equations governing fluid motion (Equation (2)) [17]. \[\frac{\partial}{\partial X_{i}}U_{i}A_{i}=0 \tag{1}\]  
\[\frac{\partial U_{i}}{\partial t}+\frac{1}{V_{f}}\left(U_{i}A_{j}\frac{ \partial u}{\partial X_{j}}\right)=-\frac{1}{\rho}\frac{\partial P}{\partial X _{i}}+G_{i}+f_{i} \tag{2}\]  
\[\rho V_{f}f_{i}=\tau_{b,i}-\left[\frac{\partial}{\partial X_{j}}\left(A_{j}S_ {ij}\right)\right];S_{ii}=-2\mu_{tot}\left[\frac{\partial U_{i}}{\partial X_{ j}}\right];S_{ij}=-\mu_{tot}\left[\frac{\partial U_{i}}{\partial X_{j}}+\frac{ \partial U_{j}}{\partial X_{i}}\right] \tag{3}\]  
where \(P\) indicates pressure, \(U_{i}\) is the mean velocity, time is denoted by \(t\), \(A_{i}\) is the fractional open area open to flow in the i direction, \(V_{f}\) is the fractional volume open to flow, \(G_{i}\) represents the body accelerations, \(f_{i}\) indicates the viscous accelerations, \(S_{ij}\) is the strain rate tensor, \(\tau_{b,i}\) denotes the wall shear stress, \(\rho\) is the density of water, \(\mu_{tot}\) is the total dynamic viscosity, which includes the effects of turbulence (\(\mu_{tot}=\mu+\mu_{T}\)), \(\mu\) is the dynamic viscosity, and \(\mu_{T}\) is the eddy viscosity. Sediment transport modeling requires accurate estimates of the near-wall shear stresses, so a good turbulence model should be selected for turbulence flows. In FLOW-3D, there are six turbulence models available: the Prandtl mixing length model, the one-equation, the two-equation k-\(\epsilon\), RNG, and k-\(\omega\) models, and a large eddy simulation, LES, model. The Renormalization Group (RNG) turbulence model [18] is used in this study. The RNG model uses equations similar to the equations for the k-\(\epsilon\) model. However, equation constants that are found empirically in the standard k-\(\epsilon\) model are derived explicitly in the RNG model. In particular, the RNG model is known to describe low-intensity turbulence flows and flows having strong shear regions more accurately. RNG turbulence model equations are shown below. \[\frac{\partial(\rho k)}{\partial t}+\frac{\partial(\rho ku_{i})}{\partial x _{i}}=\frac{\partial}{\partial x_{j}}\left[\left(\mu+\frac{\mu_{t}}{\sigma_{ k}}\right)\frac{\partial k}{\partial x_{j}}\right]+P_{k}-\rho\epsilon \tag{4}\]  
\[\frac{\partial(\rho\epsilon)}{\partial t}+\frac{\partial(\rho\epsilon u_{i})}{ \partial x_{i}}=\frac{\partial}{\partial x_{j}}\left[\left(\mu+\frac{\mu_{t}} {\sigma_{k}}\right)\frac{\partial(\epsilon)}{\partial x_{j}}\right]+C_{1} \epsilon\frac{\epsilon}{k}P_{k}-C_{2\epsilon}\rho\frac{\epsilon^{2}}{k} \tag{5}\]where \(k\) is the turbulence kinetic energy (henceforth, TKE), \(\varepsilon\) is the dissipation rate, \(\rho\) is density, \(t\) is time, \(x_{i}\) is the coordinate in the (i) axis, \(\mu\) is the dynamic viscosity, \(\mu_{t}\) is the turbulent dynamic viscosity, and \(P_{k}\) is the production of TKE. The default values of \(C_{1\varepsilon}\), \(C_{2\varepsilon}\), and \(\sigma_{k}=\sigma_{\varepsilon}\) are 1.42, 1.68, and 1.39, respectively. The turbulent viscosity is computed using the parameter \(C_{\mu}=0.85\) using Equation (6). \[\mu_{\rm T}=\rho C_{\mu}\frac{k^{2}}{\varepsilon} \tag{6}\]  
Sediment is entrained by picking-up and resuspension due to shearing and small eddies at the packed sediment interface.