## 2 Methodology  
### Hydrodynamic Model  
This study resolves 3D flow structures in meandering channels using the 3D RANS equations, which are solved with OpenFOAM, an open-source computational fluid dynamics library based on the finite volume method. The OpenFOAM-based RANS models have been extensively used for analysis of open channel flows including bed erosion [28; 38], hydraulic jump [39; 40; 41], spillway jet [41], flows around hydraulic structures [42; 43; 44], and coupled flows with porous media [45]. For incompressible fluid flows, the 3D RANS model adopts the Reynold approximation of continuity and momentum equations as follows:  
\[\frac{\partial U_{i}}{\partial x_{i}}=0 \tag{1}\]  
\[\frac{\partial U_{i}}{\partial t}+U_{j}\frac{\partial U_{i}}{\partial x_{j}}= -\frac{1}{\rho}\frac{\partial P}{\partial x_{i}}+\frac{\partial}{\partial x_ {j}}\bigg{[}v\bigg{(}\frac{\partial U_{i}}{\partial x_{j}}+\frac{\partial U_{ j}}{\partial x_{i}}\bigg{)}\bigg{]}-\frac{\partial}{\partial x_{j}}\overline{u_{i}u_{j}} \tag{2}\]  
where \(U_{i}\) is the time-averaged velocity in the Cartesian coordinate; \(P\) is the pressure; \(\rho\) is the fluid density; \(\
u\) is the kinematic viscosity; \(u_{i}\) is the velocity fluctuation; and \(\overline{u_{i}u_{j}}\) is the Reynolds stress. Under the Boussinesq approximation, \(\overline{u_{i}u_{j}}\) is defined as:  
\[-\overline{u_{i}u_{j}}=\
u_{t}\bigg{(}\frac{\partial U_{i}}{\partial x_{j}}+ \frac{\partial U_{j}}{\partial x_{i}}\bigg{)}+\frac{2}{3}k\delta_{ij} \tag{3}\]  
where \(\
u_{t}\) is the turbulent eddy viscosity; \(k\) is the turbulent kinetic energy; \(\delta_{ij}\) is the Kronecker delta. To close Equation (2) by resolving \(\
u_{t}\), this study employs two-equation turbulence closure models such as \(k-\varepsilon\), \(k-\omega\) and SST \(k-\omega\) models available with the Open FOAM library. In the \(k-\varepsilon\) model, \(\
u_{t}\) is modeled with the turbulent kinetic energy, \(k\) and the turbulent energy dissipation, \(\varepsilon\) using the relation described as:  
\[\
u_{t}=\rho C_{\mu}\frac{k^{2}}{\varepsilon} \tag{4}\]  
where \(C_{\mu}\) is the empirical constant. The additional transport equations for simulating \(k\) and \(\varepsilon\) are given as:  
\[\frac{\partial k}{\partial t}+\frac{\partial(U_{i}k)}{\partial x_{i}}=\frac {\partial}{\partial x_{j}}\bigg{[}\bigg{(}\frac{\
u_{t}}{\sigma_{k}}+v\bigg{)} \frac{\partial k}{\partial x_{j}}\bigg{]}-\varepsilon+G_{k} \tag{5}\]  
\[\frac{\partial\varepsilon}{\partial t}+\frac{\partial(U_{i}\varepsilon)}{ \partial x_{i}}=\frac{\partial}{\partial x_{j}}\bigg{[}\bigg{(}\frac{\
u_{t }}{\sigma_{\varepsilon}}+v\bigg{)}\frac{\partial\varepsilon}{\partial x_{j} }\bigg{]}+C_{\varepsilon 1}G_{k}\frac{\varepsilon}{k}-C_{\varepsilon 2}\frac{ \varepsilon^{2}}{k} \tag{6}\]  
where \(\sigma_{k}\), \(\sigma_{\varepsilon}\), \(C_{\varepsilon 1}\), and \(C_{\varepsilon 2}\) are the empirical constants; and \(G_{k}\) is the production of turbulent kinetic energy due to the mean velocity gradient. The empirical constants of the \(k-\varepsilon\) model have been reported as \(C_{\mu}=0.09\), \(C_{\varepsilon 1}=1.44\), \(C_{\varepsilon 2}=1.92\), \(\sigma_{k}=1.0\), \(\sigma_{\varepsilon}=1.3\) at high Reynolds number [46]. The \(k-\varepsilon\) model is well known as the simplest turbulence model providing acceptable accuracy for diverse hydraulic problems with good convergence. However, the inherent limitations originated from the turbulent viscosity relation in Equation (4) and the \(\varepsilon\) equation lead to inaccurate simulations for flows near the boundary layers. Therefore, the damping function is required to correct the overestimated \(\
u_{t}\) within the boundary flow [47].