*3.1. 2D Numerical Model*

The conservation form of the governing equations is expressed as follows:

$$\frac{\partial \mathbf{Q}\_{\rm v}}{\partial t} + \frac{\partial \mathbf{F}\_{\rm adj}}{\partial \mathbf{x}} + \frac{\partial \mathbf{G}\_{\rm adj}}{\partial y} = \frac{\partial \mathbf{F}\_{\rm diff}}{\partial \mathbf{x}} + \frac{\partial \mathbf{G}\_{\rm diff}}{\partial y} + \mathbf{S}\_{\rm s} \tag{1}$$

Here,

$$\begin{aligned} \mathbf{Q}\_{\mathbf{v}} &= \begin{bmatrix} h \\ h\mathbf{u} \\ h\mathbf{v} \\ h\mathbf{C} \end{bmatrix}, \mathbf{F}\_{\text{adv}} = \begin{bmatrix} h\mathbf{u} \\ h\mathbf{u}^2 + gh^2/2 \\ h\mathbf{u}\mathbf{v} \\ h\mathbf{u}\mathbf{C} \end{bmatrix}, \mathbf{G}\_{\text{adv}} = \begin{bmatrix} h\mathbf{v} \\ h\mathbf{u}\mathbf{v} \\ h\mathbf{v}^2 + gh^2/2 \\ h\mathbf{v}\mathbf{C} \end{bmatrix}, \\ \mathbf{F}\_{\text{diff}} &= \begin{bmatrix} 0 \\ hT\_{\text{xx}}/\rho \\ hT\_{\text{xy}}/\rho \\ D\_x\text{\ddot{\eta}}(h\mathbf{C})/\beta \mathbf{x} \end{bmatrix}, \mathbf{G}\_{\text{diff}} = \begin{bmatrix} 0 \\ hT\_{\text{xy}}/\rho \\ hT\_{\text{yy}}/\rho \\ D\_y\text{\ddot{\eta}}(h\mathbf{C})/\beta y \end{bmatrix}, \mathbf{S}\_{\text{s}} = \begin{bmatrix} 0 \\ gh(\mathbf{s}\_{\text{0x}} - \mathbf{s}\_{\text{f}\mathbf{x}}) \\ gh(\mathbf{s}\_{\text{0y}} - \mathbf{s}\_{\text{f}\mathbf{y}}) \\ q\mathbf{s} - q\_{\text{sd}} \end{bmatrix} \end{aligned} (2)$$

where Qv is the conserved physical vector; Fadv and Gadv are the advection flux vectors in the *x* and *y* directions, respectively; Fdiff and Gdiff are the diffusion flux vectors in the *x* and *y* directions, respectively; S*s* is the source term; *h* is the water depth; *u* and *v* are the depth-averaged velocity components in the *x* and *y* directions, respectively; ρ is the density of water; *Txx*, *Txy*, and *Tyy* are the depth-averaged turbulent stresses; *C* is the depth-averaged suspended sediment concentration; *g* is the gravitational acceleration; D*x* and <sup>D</sup>*y* are the diffusion coefficients. *s*0*x* and *<sup>s</sup>*0*y* are the bed slopes in the *x* and *y* directions, respectively, and they are expressed as:

$$\mathbf{s}\_{0\mathbf{x}} = -\partial \mathbf{z}\_{b} / \partial \mathbf{x} \text{ and } \mathbf{s}\_{0y} = -\partial \mathbf{z}\_{b} / \partial y \tag{3}$$

*zb* is the bed elevation; *s*f*x* and *<sup>s</sup>*f*y* are the friction slopes in the *x* and *y* directions, respectively; and *Txx*, *Txy*, and *Tyy* are the depth-averaged turbulent stresses. On the basis of Manning formulas, the friction slopes can be expressed as:

$$s\_{\rm fx} = \frac{\imath m\_{\rm M}^2 \sqrt{\imath^2 + \upsilon^2}}{h^{\frac{4}{3}}}, \; s\_{\rm fy} = \frac{\imath m\_{\rm M}^2 \sqrt{\imath^2 + \upsilon^2}}{h^{\frac{4}{3}}} \tag{4}$$

where *n*M is the Manning roughness coefficient.

According to the particle size ranges, suspended sediment particles can be divided into cohesive and non-cohesive particles [21]. For cohesive sediment, the deposition flux (depositional rate) can be determined using the following relation [22].

$$q\_{sd} = P w\_s \text{C\\_P = 1} - \sqrt{\tau\_{bx}^2 + \tau\_{ly}^2} / \tau\_{cd} \tag{5}$$

where *P* is the probability of deposition and τ*cd* is the critical shear stress of deposition. Moreover, *ws* is the settling velocity and defined as:

$$w\_s = 250 \, d\_m^{0.2} \tag{6}$$

where *dm* is the median particle size (μm). The entrainment flux (erosion rate) for cohesive materials was given by Liu et al. (2002) [23] and can be computed from:

$$q\_{\rm sc} = \mathbb{E}(\sqrt{\tau\_{bx}^2 + \tau\_{by}^2}/\tau\_{cx} - 1) \tag{7}$$

where *E* is the erosion parameter and τ*ce* is the critical shear stress of erosion. In general, the critical shear stresses of deposition and erosion differ in various study areas. Based on the experimental data and model calibration, Liu et al. (2002) [23] provided the suggested values of <sup>τ</sup>*cd*(= 0.05 <sup>N</sup>/m2) and <sup>τ</sup>*ce*(<sup>=</sup> 0.1 <sup>N</sup>/m2) for modeling cohesive sediment transport in Equations (5) and (7), respectively.

For the numerical discretization in the framework of the finite volume method, the finite-volume multi-stage (FMUSTA) scheme was adopted to estimate the numerical flux through each cell interface in order to solve Equation (1) [24]. In the treatment of source terms, the hydrostatic reconstruction method was employed to prevent the imbalance between flux gradients and source terms. The presented 2D numerical model is practically suitable for complex flow conditions (mixed regimes) and is appropriate for special problems related to wetting and drying. A detailed description of the numerical formulation was provided by Guo et al. (2011) [18].
