2.1.1. Governing Equations

The governing equations are the two-dimensional incompressible Navier-Stokes, that are composed by the momentum conservation:

$$\frac{\partial u}{\partial t} + \frac{\partial uu}{\partial \mathbf{x}} + \frac{\partial uv}{\partial z} = \quad -\frac{\partial p}{\partial \mathbf{x}} + \nu \left( \frac{\partial^2 u}{\partial \mathbf{x}^2} + \frac{\partial^2 u}{\partial z^2} \right),$$

$$\frac{\partial w}{\partial t} + \frac{\partial wu}{\partial \mathbf{x}} + \frac{\partial ww}{\partial z} = \quad -\frac{\partial p}{\partial z} + \nu \left( \frac{\partial^2 w}{\partial \mathbf{x}^2} + \frac{\partial^2 w}{\partial z^2} \right) - \mathbf{g},\tag{1}$$

and the incompressibility condition:

$$\frac{\partial u}{\partial x} + \frac{\partial w}{\partial z} = 0,\tag{2}$$

where *u*, *w* are the components of the velocity field in the *x*− and *z*− direction, *ν* = *μ*/*ρ* is the kinematic viscosity coefficient, *μ* is the dynamic viscosity coefficient, *ρ* is the fluid density, *p* is the pressure term normalized by the fluid density, and *g* is the gravity acceleration. If the pressure is assumed to be locally hydrostatic, it can be expressed in terms of the atmospheric pressure *pa* and the hydraulic head *η* as:

$$p = p\_a + \lg(\eta - z). \tag{3}$$

By means of Equation (3) it is possible to rewrite the momentum Equations (1) in terms of the hydraulic head *η*:

$$\frac{\partial u}{\partial t} + \frac{\partial uu}{\partial x} + \frac{\partial uw}{\partial z} = -g\frac{\partial \eta}{\partial x} + \nu \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial z^2}\right),$$

$$\frac{\partial w}{\partial t} + \frac{\partial wu}{\partial x} + \frac{\partial ww}{\partial z} = -g\frac{\partial \eta}{\partial z} + \nu \left(\frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial z^2}\right). \tag{4}$$

Further introducing the mass conservation law for the suspended sediment concentration *c* (here defined as volume concentration, hence simplifying the mass conservation law assuming constant sediment density), and defining its relation with the sediment level *S* = *S*(*x*, *z*, *t*), two more equations need to be added to the system:

$$\frac{\partial \mathbf{c}}{\partial t} + \frac{\partial \mathbf{u} \mathbf{c}}{\partial \mathbf{x}} + \frac{\partial (\mathbf{w} - \mathbf{w}\_c) \mathbf{c}}{\partial z} = f\_{c\mathbf{S}}(\mathbf{S}, \mathbf{c}) \tag{5}$$

$$\frac{\partial S}{\partial t} = f\_{\text{Sc}}(S, c), \tag{6}$$

where *wc* is the falling velocity of the suspended sediment; *fSc* and *fcS* are two functions that regulate the mass transfer between the passive concentration *c* and the sediment level *S*. Note that in Equations (5) and (6) the sediment dynamics are intrinsically assumed to be driven only by the sedimentation/erosion process through the functions *fSc* and *fcS*. The suspended concentration is expressed as volume of suspended sediment per unit of control volume (i.e., volume concentration relative to each computational cell), thus it can range in *c* ∈ [0, 1 − *φ*], where *φ* indicates the sediment porosity. Similarly, the sediment level *S* is defined per unit of volume, thus it can vary in the range *S* ∈ [0, 1], where *S* = 1 corresponds to a computational cell completely filled by the sediment. Note that the suspended sediment is assumed to be passively transported by the fluid, while the deposited sediment level actively interacts with the fluid dynamics as in general it may vary in time, thus changing the boundaries of the fluid domain. The exchange between suspended and deposited sediment must satisfy the mass conservation law:

$$\frac{d}{dt}\left[c + S(1 - \phi)\right] = 0.\tag{7}$$

.

## 2.1.2. Staggered Mesh

The Partial Differential Equations (PDE) system composed by Equations (2), (4) and (5) is numerically solved in a domain <sup>Ω</sup>(*t*) ⊆ <sup>Ω</sup> ∈ <sup>2</sup> discretized by a regular mesh <sup>Ω</sup> = *i*,*k* Ω*i*,*k*, where <sup>Ω</sup>*i*,*<sup>k</sup>* = [*xi*<sup>−</sup> <sup>1</sup> 2 , *xi*<sup>+</sup> <sup>1</sup> 2 ] <sup>×</sup> [*zk*<sup>−</sup> <sup>1</sup> 2 , *zk*<sup>+</sup> <sup>1</sup> 2 ] are rectangular control volumes centered in (*xi*, *zk*). The mesh is assumed to be homogeneous with Δ*x* = *xi*<sup>+</sup> <sup>1</sup> 2 <sup>−</sup> *xi*<sup>−</sup> <sup>1</sup> 2 and Δ*z* = *zk*<sup>+</sup> <sup>1</sup> 2 <sup>−</sup> *zk*<sup>−</sup> <sup>1</sup> 2 , for *i* = 1, ... , *Imax* and *k* = 1, ... , *Kmax*. The main quantities are defined following a staggered approach. In particular, the hydraulic head *η*, the suspended concentration *c* and the sediment level *S* are defined in the cell centers, namely *η* = *ηik*, *c* = *cik*, *S* = *Sik*. The horizontal velocity is defined in the center of the vertical edges, while the vertical velocity in the center of horizontal edges, namely, *<sup>u</sup>* <sup>=</sup> *ui*<sup>±</sup> <sup>1</sup> <sup>2</sup> ,*<sup>k</sup>* and *<sup>w</sup>* <sup>=</sup> *wi*,*k*<sup>±</sup> <sup>1</sup>

2 The system is solved for discrete times *t* = *t* 0, *t* 1, ... , *t <sup>n</sup>*, ..., hence the solution at a certain time *t n* is marked with the upper index *n*, namely *ξ*(*t <sup>n</sup>*) := *ξn*, while the time interval between two consecutive steps is Δ*t <sup>n</sup>* = *t <sup>n</sup>*+<sup>1</sup> <sup>−</sup> *<sup>t</sup> <sup>n</sup>*. Note that the effective volume in each control cell can have different values due to the presence of the sediment, and can be computed as *Vi*,*<sup>k</sup>* = Δ*x*Δ*z*(1 − *Sik*). Since *Si*,*<sup>k</sup>* ∈ [0, 1], the volume computed in this way satisfies 0 ≤ *Vi*,*<sup>k</sup>* ≤ Δ*x*Δ*z*. The effective edge length can change as well, being defined above the sediment top: for any time *t <sup>n</sup>* the effective vertical and horizontal edges lengths are defined as Δ*z<sup>n</sup> <sup>i</sup>*<sup>±</sup> <sup>1</sup> <sup>2</sup> ,*<sup>k</sup>* and <sup>Δ</sup>*x<sup>n</sup> <sup>i</sup>*,*k*<sup>±</sup> <sup>1</sup> 2 , see Figure 1. For non-saturated cells (i.e., *S* < 1), the effective edges lengths are computed using an upwind procedure:

$$
\Delta z\_{i+\frac{1}{2},k}^{n} = \begin{cases}
\frac{V\_{i,k}^{n}}{\Delta x} & u\_{i+\frac{1}{2},k}^{n} > 0 \\
\frac{V\_{i+1,k}^{n}}{\Delta x} & u\_{i+\frac{1}{2},k}^{n} < 0 \\
\frac{\max\{V\_{i,k'}^{n}, V\_{i+1,k}^{n}\}}{\Delta x} & u\_{i+\frac{1}{2},k}^{n} = 0
\end{cases}
\quad \Delta x\_{i,k+\frac{1}{2}}^{n} = \begin{cases} 0 & S\_{i,k} = 1 \\\\ \Delta x & S\_{i,k} \neq 1 \end{cases} \tag{8}
$$

We further assume that Δ*z<sup>n</sup> <sup>i</sup>*<sup>±</sup> <sup>1</sup> <sup>2</sup> ,*<sup>k</sup>* <sup>=</sup> <sup>Δ</sup>*x<sup>n</sup> <sup>i</sup>*,*k*<sup>±</sup> <sup>1</sup> 2 = 0 for any (*i*, *k*) where *Si*,*<sup>k</sup>* = 1.

**Figure 1.** Schematic of the used staggered mesh.
