2.1.3. Numerical Method

In order to avoid the development of a fully nonlinear system, the convective and the viscous terms in Equation (4) are discretized explicitly, while the velocity field and the hydraulic head in Equations (2) and (4) are discretized implicitly. A finite volume approximation of the continuity Equation (2) reads:

$$
\Delta z\_{i + \frac{1}{2}k}^n u\_{i + \frac{1}{2}k}^{n + 1} - \Delta z\_{i - \frac{1}{2}k}^n u\_{i - \frac{1}{2}k}^{n + 1} + \Delta x\_{i, k + \frac{1}{2}}^n w\_{i, k + \frac{1}{2}}^{n + 1} - \Delta x\_{i, k - \frac{1}{2}}^n w\_{i, k - \frac{1}{2}}^{n + 1} = 0. \tag{9}$$

Consistently, a finite difference approximation for the momentum Equations (4) reads:

$$u\_{i+\frac{1}{2},k}^{n+1} = F u\_{i+\frac{1}{2},k}^n - \Delta t \lg \frac{\eta\_{i+1,k}^{n+1} - \eta\_{i,k}^{n+1}}{\Delta x} \tag{10}$$

$$w\_{i,k+\frac{1}{2}}^{n+1} = \mathcal{F}w\_{i,k+\frac{1}{2}}^n - \Delta t \lg \frac{\eta\_{i,k+1}^{n+1} - \eta\_{i,k}^{n+1}}{\Delta z},\tag{11}$$

where *Fu<sup>n</sup> i*+ <sup>1</sup> 2 ,*k* , *Fw<sup>n</sup> i*,*k*+ <sup>1</sup> 2 contain all the explicit contributions of the convective and viscous terms, that in this paper are expressed through the explicit conservative formulation proposed in [39]. Substituting the discrete velocities *un*+<sup>1</sup> *<sup>i</sup>*<sup>±</sup> <sup>1</sup> 2 ,*k* , *wn*+<sup>1</sup> *<sup>i</sup>*,*k*<sup>±</sup> <sup>1</sup> 2 as expressed in Equations (10) and (11) into Equation (9), we obtain a linear system for the unknown hydraulic head:

$$\frac{\Delta t \mathcal{g}}{\Delta x} \left[ \Delta z\_{i + \frac{1}{2}, k}^{n} \left( \eta\_{i + 1, k}^{n+1} - \eta\_{i, k}^{n+1} \right) - \Delta z\_{i - \frac{1}{2}, k}^{n} \left( \eta\_{i, k}^{n+1} - \eta\_{i - 1, k}^{n+1} \right) \right]$$

$$+ \frac{\Delta t \mathcal{g}}{\Delta z} \left[ \Delta x\_{i, k + \frac{1}{2}}^{n} \left( \eta\_{i, k + 1}^{n+1} - \eta\_{i, k}^{n+1} \right) - \Delta x\_{i, k - \frac{1}{2}}^{n} \left( \eta\_{i, k}^{n+1} - \eta\_{i, k - 1}^{n+1} \right) \right] = b\_{i, k}^{n} \tag{12}$$

where the term *b<sup>n</sup> <sup>i</sup>*,*<sup>k</sup>* contains all the known terms evaluated at the time *t n*:

$$b\_{i,k}^{n} = \Delta z\_{i + \frac{1}{2},k}^{n} F u\_{i + \frac{1}{2},k}^{n} - \Delta z\_{i - \frac{1}{2},k}^{n} F u\_{i - \frac{1}{2},k}^{n} + \Delta x\_{i,k + \frac{1}{2}}^{n} F w\_{i,k + \frac{1}{2}}^{n} - \Delta x\_{i,k - \frac{1}{2}}^{n} F w\_{i,k - \frac{1}{2}}^{n} \tag{13}$$

Equation (12) is a five-points diagonal system, symmetric and at least semi-positive definite (e.g., [36]). Therefore, it can be solved using a matrix-free implementation of the conjugate gradient method. Once the hydraulic head *ηn*+<sup>1</sup> is known, the velocity field can be easily updated through the explicit formulae in Equations (10) and (11).

As for the scalar transport, similarly to [38], we refer to a semi-implicit finite volume approximation based on upwind fluxes:

$$c\_{i,k}^{\*} = c^{n} + \frac{\Delta t}{\Delta x \Delta z} \left[ \Delta z\_{i + \frac{1}{2}k}^{n} c\_{i + \frac{1}{2}k}^{n, \mu p} \left( u\_{i + \frac{1}{2}k}^{n+1} \right) - \Delta z\_{i - \frac{1}{2}k}^{n} c\_{i - \frac{1}{2}k}^{n, \mu p} \left( u\_{i - \frac{1}{2}k}^{n+1} \right) \right]$$

$$\Delta x\_{i, k + \frac{1}{2}}^{n} c\_{i, k + \frac{1}{2}}^{n, \mu p} \left( w\_{i, k + \frac{1}{2}}^{n + 1} - w\_{\delta} \right) - \Delta x\_{i, k - \frac{1}{2}}^{n} c\_{i, k - \frac{1}{2}}^{n, \mu p} \left( w\_{i, k - \frac{1}{2}}^{n + 1} - w\_{\delta} \right) \right],\tag{14}$$

where *c n*,*up* · is the upwind contribution defined as:

$$c\_{i+\frac{1}{2},k}^{n,up}(V) = \frac{1}{2} \left[ c\_{i,k}^n(|V|+V) + c\_{i+1,k}^n(|V|-V) \right] \tag{15}$$

$$c\_{i,k+\frac{1}{2}}^{n,\mu p}(V) = \frac{1}{2} \left[ c\_{i,k}^{n}(|V|+V) + c\_{i,k+1}^{n}(|V|-V) \right] \tag{16}$$

If there is no sedimentation/erosion, then *cn*+<sup>1</sup> *<sup>i</sup>*,*<sup>k</sup>* = *c*<sup>∗</sup> *<sup>i</sup>*,*k*, otherwise the mass transfer is performed assuming that the fluid is in a local-static equilibrium. First, the dynamic part of the suspended sediment is computed through Equation (14), then the result is used to update explicitly the new sediment level through the discrete version of Equation (6):

$$S\_{i,k}^{n+1} = S\_{i,k}^{n} + \frac{\Delta t}{\Delta x \Delta z} Q\_{SC}(S\_{i,k'}^{n} c\_{i,k}^{\*}), \qquad \text{where} \qquad Q\_{SC}(S\_{i,k'}^{n} c\_{i,k}^{\*}) = \int\_{\Omega\_{i,k}} f\_{SC}(S\_{i,k'}^{n} c\_{i,k}^{\*}).\tag{17}$$

The sediment level at the time step *n* + 1 as resulting from Equation (17) allows for updating the volumes/edges lengths specified in Equation (8). Finally, the concentration at the time step *n* + 1 is updated according to the discrete version of the mass conservation law (7) that reads:

$$c\_{i,k}^{n+1} = c\_{i,k}^{\*} + \left(S\_{i,k}^{n} - S\_{i,k}^{n+1}\right)(1 - \phi). \tag{18}$$

We note that by substituting Equation (14) into (18) and integrating *cn*+<sup>1</sup> *<sup>i</sup>*,*<sup>k</sup>* over the water column, one obtains the discrete mass conservation law for the suspended transport, where the source term is the discrete version of *<sup>∂</sup>zb <sup>∂</sup><sup>t</sup>* (<sup>1</sup> <sup>−</sup> *<sup>φ</sup>*), with *zb* total sediment level. However, the form in (14)–(18) is more general since it does not require a single value function *zb* = *zb*(*x*, *t*), but it can possibly be applied to cases with multiple sediment/water interfaces.

For the chosen explicit discretization of the nonlinear convective terms, the method is stable under a Courant–Friedrichs–Lewy (CFL) type restriction based on the fluid velocity (e.g., [40]),

$$
\Delta t \le \frac{1}{\frac{|u|}{\Delta x} + \frac{|w|}{\Delta z} + \frac{2\nu}{\Delta x^2} + \frac{2\nu}{\Delta z^2}}.\tag{19}
$$

The method becomes unconditionally stable if an Eulerian-Lagrangian scheme is adopted, in combination with a sub-step approach for the evolution of the concentration and sediment level [37,41].
