2.2.2. Numerical Discretizations

The adopted HDM uses a θ semi-implicit formulation [24–26], while finite-volume and finite-difference methods are combined. Momentum equations are solved within a finite-difference framework and using operator-splitting techniques. The θ semi-implicit method is used to advance the time stepping. Correspondingly, the gradient of the free-surface elevation is discretized into explicit and implicit parts. A point-wise Eulerian-Lagrangian method (ELM), using the multistep backward

Euler technique [23,27], is used to solve the advection term. The horizontal diffusion term is discretized using an explicit center-difference method.

When the advection term is solved by the ELM, the velocities are updated at once and are denoted by *ubt* and *vbt*. The horizontal momentum equations in the local horizontal *x*-, *y*-directions of unstructured grids are then discretized as follows (at side *j*)

$$\left(1+\Delta t \text{g} \boldsymbol{u}\_{\text{m}}^{2} \frac{\sqrt{\boldsymbol{\mu}\_{\text{it},j}^{n}+\boldsymbol{v}\_{\text{it},j}^{n}}}{\boldsymbol{h}\_{j}^{n4/3}}\right) \boldsymbol{u}^{n+1} = \boldsymbol{u}\_{\text{it},j}^{n} - \Delta t \text{g} \Big[ (1-\theta) \frac{\eta\_{i(j,2)}^{n} - \eta\_{i(j,1)}^{n}}{\delta\_{j}} + \theta \frac{\eta\_{i(j,2)}^{n+1} - \eta\_{i(j,1)}^{n+1}}{\delta\_{j}} \Big] + \Delta t \text{E} \boldsymbol{x}\_{j}^{n} \tag{7a}$$

$$\left(1+\Delta t \text{g} \boldsymbol{u}\_{m}^{2} \frac{\sqrt{\boldsymbol{u}\_{\mathrm{b},j}^{n}+\boldsymbol{v}\_{\mathrm{b},j}^{n}}}{\boldsymbol{h}\_{j}^{n+1/3}}\right) \mathbf{v}^{n+1} = \mathbf{v}\_{\mathrm{b},j}^{n} - \Delta t \text{g} \left[ (1-\theta) \frac{\eta\_{\mathrm{ip}(j,2)}^{n} - \eta\_{\mathrm{ip}(j,1)}^{n}}{\boldsymbol{L}\_{j}} + \theta \frac{\eta\_{\mathrm{ip}(j,2)}^{n+1} - \eta\_{\mathrm{ip}(j,1)}^{n+1}}{\boldsymbol{L}\_{j}} \right] + \Delta t E\_{Y} \boldsymbol{\upmu}^{n} \tag{7b}$$

where θ is the implicit factor and Δ*t* the time step; superscripts "*n*" indicate the *n*-th time step; for simplicity, the explicitly discretized horizontal diffusion term is not expanded here; the riverbed friction is discretized using *ubt* and *vbt* to enhance computation stability. The explicitly discretized horizontal diffusion term is not expanded here for simplicity, and denoted by *EX* and *EY* in x- and y-directions, respectively. The η at nodes is regarded as auxiliary variables, which are interpolated from water-level values of neighboring cells.

When explicit terms of the discretized momentum equations are incorporated, the unknowns (free-surface elevation η) emerge. Equation (7a,b) are then transformed into (at side *j*)

$$\mu\_{j}^{n+1} = G\_{j}^{n} / A\_{j}^{n} - \Theta \text{g} \Delta t \frac{\eta\_{i(j,2)}^{n+1} - \eta\_{i(j,1)}^{n+1}}{\delta\_{j}} / A\_{j}^{n} \tag{8a}$$

$$
\sigma\_{j}^{n+1} = F\_{j}^{n} / A\_{j}^{n} - \theta \lg \frac{\eta\_{ip(j,2)}^{n+1} - \eta\_{ip(j,1)}^{n+1}}{L\_{j}} / A\_{j}^{n} \tag{8b}
$$

where *Anj* = 1 + <sup>Δ</sup>*tgn*<sup>2</sup>*m ubtnj* 2 + *ubtnj* 2/*hnj* 4/3; *Gnj* , *Fnj* are the incorporated explicit terms respectively in the horizontal *x*-, *y*-directions, *Gnj* = *<sup>u</sup>nbt*,*<sup>j</sup>* − Δ*tg*(1 − θ) <sup>η</sup>*ni*(*j*,<sup>2</sup>)<sup>−</sup><sup>η</sup>*ni*(*j*,<sup>1</sup>) δ*j* + Δ*tEXnj* , *Fnj* = *<sup>v</sup>nbt*,*<sup>j</sup>* − Δ*tg*(1 − θ) <sup>η</sup>*nip*(*j*,<sup>2</sup>)<sup>−</sup><sup>η</sup>*nip*(*j*,<sup>1</sup>) *Lj* + Δ*tEYnj* .

To achieve good mass conservation, the depth-integrated continuity equation, Equation (1), is discretized by the finite-volume method, which is given by (at cell *i*)

$$P\_i \eta\_i^{n+1} = P\_i \eta\_i^n - \Theta \Delta t \sum\_{l=1}^{34(i)} s\_{iJ} L\_{j(iJ)} h\_{j(iJ)}^n u\_{j(iJ)}^{n+1} - (1 - \Theta) \Delta t \sum\_{l=1}^{34(i)} s\_{iJ} L\_{j(iJ)} h\_{j(iJ)}^n u\_{j(iJ)}^n \tag{9}$$

where *l* is the side index of cell *i*, and *l* = 1, 2, ... , *i34*(*i*).

The velocity–pressure coupling is performed by substituting *uj<sup>n</sup>*+<sup>1</sup> and *vj<sup>n</sup>*+<sup>1</sup> of Equation (8a,b) into the discrete depth-integrated continuity equation. This substitution results in a wave propagation equation with cell water levels (η) as unknowns. Using the topology relations among the cells, the resulting discrete wave propagation equation is given by (at cell *i*)

$$\begin{aligned} &\quad P\_{i}\eta\_{i}^{n+1} + \mathcal{g}\theta^{2}\Delta t^{2} \sum\_{l=1}^{i34(i)} \frac{L\_{j(iJ)}}{\delta\_{\dot{j}}} \mathcal{h}\_{j(iJ)}^{n} \Big(\eta\_{i}^{n+1} - \eta\_{i i \mathcal{G}(iJ)}^{n+1}\Big) / A\_{j(iJ)}^{n} \\ &= P\_{i}\eta\_{i}^{n} - \theta\Delta t \sum\_{l=1}^{i34(i)} \mathbf{s}\_{l,l}L\_{j(iJ)}\mathcal{h}\_{j(iJ)}^{n} \mathcal{G}\_{j(iJ)}^{n} / A\_{j(iJ)}^{n} - (1 - \theta)\Delta t \sum\_{l=1}^{i} \mathbf{s}\_{l,l}L\_{j(iJ)}\mathcal{h}\_{j(iJ)}^{n} d\_{j(iJ)}^{n} \end{aligned} \tag{10}$$

The HDM solves the vertically averaged 2D shallow water equations at three steps. First, all the explicit terms (advection, diffusion, riverbed friction, and the explicit part of free-surface

gradients) in momentum equations are explicitly computed to obtain the provisional velocities. Second, the velocity-pressure coupling is performed by substituting the expressions of normal velocity components into the discrete continuity equation, where a wave propagation equation is constructed and solved to obtain new water levels. Third, a back substitution of the new water levels into the momentum equations is performed to ge<sup>t</sup> the final velocity field.

For each fraction of nonuniform sediment, one transport equation must be solved. The STM is advanced fully explicitly, and the transport equation is discretized as (for fraction *k*)

$$\mathbf{C}\_{\mathbf{k},i}^{n+1} = \mathbf{C}\_{\mathbf{k},\mathbf{l};i}^{n+1} + \frac{\Delta t}{P\_{\mathbf{l}}\mathbf{l}\_{\mathbf{i}}^{n}} \sum\_{l=1}^{\text{D4}(i)} \left\{ \mathbf{s}\_{i\mathbf{l}} L\_{j(i\mathbf{l})} \mathbf{h}\_{j(i\mathbf{l})}^{n} \left[ \left( \frac{\nu\_{\text{I}}}{\sigma\_{\text{c}}} \right)\_{j(i\mathbf{l})}^{\text{H}} \frac{\mathbf{C}\_{\mathbf{k},\mathbf{i}[j(i\mathbf{l}),2]}^{\text{H}} - \mathbf{C}\_{\mathbf{k},\mathbf{i}[j(i\mathbf{l}),1]}^{\text{H}}}{\delta\_{j(i\mathbf{l})}} \right] \right\} + \frac{\Delta t}{h\_{\text{i}}^{n}} a\_{\text{k}} \varpi\_{\text{sk}} \{ \mathbf{S}\_{\mathbf{\*}\mathbf{k}}^{n} - \mathbf{C}\_{\mathbf{k}}^{n} \} \tag{11}$$

where Δ*t* is the time step for the STM; *Cbt* is the solution to the advection subequation. The *Cbt* is calculated using a recently developed finite-volume ELM (FVELM) [28], where mass is conserved and large time steps (for which the Courant-Friedrichs-Lewy number (CFL) can be much greater than 1) are allowed.

For the FVELM, the geometrical computation which is common for each sediment fraction can be reused. When the most time-consuming parts (calculations of trajectories and interpolation weights) are avoided, only a relatively very small computation cost is added for solving each additional sediment fraction. Therefore, the FVELM allows constructing e fficient algorithms for solving the transport of a large number of sediment fractions, and this property of the FVELM is defined as the multiscalar property. Benefiting from "allowing large time steps, parallelizable, multiscalar property", the FVELM is much more e fficient than the traditional Eulerian advection schemes in solving the transport of nonuniform sediment with several fractions.
