Numerical Methods (1D FV)

Discrete Scheme for SWE. The SWE express volume and momentum conservation. The continuity equation is discretized using finite volume approximations. For the momentum equation, the type of discretization will vary depending on the term. The Crank-Nicolson method is also used to weight the contribution of variables at time steps n and n+1. However, the different nature of the equations will call for the use of a more elaborate solver scheme.

Mass Conservation

The continuity equation representing water mass (and volume) conservation is discretized as

1)	$\begin{array}{l}\displaystyle \displaystyle \frac{\Omega_i^{n+1} - \Omega_i^{n} }{\Delta t} - V_k^{n+\theta} A_k^{n+\theta} + V_{k+1}^{n+\theta} A_{k+1}^{n+\theta} = Q_i\end{array}$

where $\begin{array}{l}\Delta t\end{array}$ is the time step, $\begin{array}{l}\Omega_i\end{array}$ is the volume at cell $\begin{array}{l}i\end{array}$ , $\begin{array}{l}V_k\end{array}$ is the cross-section velocity, $\begin{array}{l}A_k\end{array}$ is the cross-sectional area, and $\begin{array}{l}Q_i\end{array}$ represents the cell sources and sinks. The cross-sectional velocity and area are computed here with $\begin{array}{l}\theta\end{array}$ -averaging. Here $\begin{array}{l}k\end{array}$ and $\begin{array}{l}k+1\end{array}$ represents the upstream and downstream cross-sections of cell $\begin{array}{l}i\end{array}$ , respectively.

Momentum Conservation

Since velocities are computed on cross-sections, the momentum equations are not located on a computational cell. The discrete equations are built based on a semi-implicit scheme in which only the acceleration, barotropic pressure gradient and bottom friction terms contain variables for which the equation is solved. Other terms of the momentum equation are still computed based on the -method, but their contribution is smaller and so they are considered explicit forcing function terms and moved to the right-hand side of the linearized system.

Acceleration. Acceleration terms are discretized using a Lagrangian approach as:

2)	$\begin{array}{l}\displaystyle \displaystyle \left ( \frac{\partial V_j}{\partial t} + V_j \frac{\partial V_j}{\partial x} \right )_k = \left ( \frac{D V_j}{D t} \right )_k = \frac{V_{k,j}^{n+1} - V_{X(j)}^n }{\Delta t}, \: \: j \in \{ch,lob,rob\}\end{array}$

where $\begin{array}{l}V_{i,j}^{n+1}\end{array}$ is the current velocity and $\begin{array}{l}V_{X(j)}^{n}\end{array}$ is evaluated at a location $\begin{array}{l}X(j)\end{array}$ . This location is found by integrating the velocity field back in time starting from the location of the computational face. Location $\begin{array}{l}X(i)\end{array}$ does not in general correspond to a cross-section, so an interpolation technique is applied.

Integration of the velocity is done in steps using the interpolated velocity field in each cell. In practice, this is equivalent to subdividing the integration time step into smaller sub-steps with a Courant number of one or less and increasing the robustness of the computation. In contrast to the explicit Eulerian framework, the semi-Lagrangian scheme allows for the use of large time steps without limiting the stability and with a much-reduced artificial diffusion (regarded as the interpolation error).

Barotropic Pressure Gradient. Recall that the momentum equation is computed at faces, but the water surface elevation term is computed at cells. This staggered grid makes the barotropic pressure gradient ideal for utilizing the simple two-point stencil described previously. In addition, the treated semi-implicitly as

3)	$\begin{array}{l}\displaystyle \displaystyle \left ( g \frac{\partial z_s}{\partial x_i} \right )_k \approx g \frac{ z_{s,i}^{n+\theta} -z_{s,i-1}^{n+\theta}}{\Delta x_{k,j}}, \: \: j \in \{ch,lob,rob\}\end{array}$

where $\begin{array}{l}g\end{array}$ is the gravitational acceleration, $\begin{array}{l}z_{s}^{n+\theta} = (1 - \theta) z_{s}^{n} + \theta z_{s}^{n+1}\end{array}$ , and $\begin{array}{l}\Delta x_{k,i}\end{array}$ is the distance between the neighboring cells along i. As a consequence of the implicit weighting, the barotropic pressure gradient consists of two parts, an implicit term with weight $\begin{array}{l}\theta\end{array}$ and an explicit term with weight $\begin{array}{l}(1-\theta)\end{array}$ .

Momentum Diffusion. Momentum diffusion represents the sum of molecular and turbulence mixing, and momentum dispersion. The momentum diffusion formulation based on the Laplacian of the current velocity is discretized as

4)	$\begin{array}{l}\displaystyle \displaystyle \left ( \nu_{t,j} \frac{\partial^2 V_j } {\partial x^2 } \right )_k \approx \nu_{t,k,j} \left ( \frac{\partial^2 V_j } {\partial x^2 } \right )_X^n, \: \: j \in \{ch,lob,rob\}\end{array}$

where is the explicit face eddy viscosity, and is the Laplacian at location X. The Laplacian is computed at nodes and spatially interpolated at the location X obtained from the acceleration advection. The Laplacian field is explicit in the numerical solution so it will depend only on values computed for the previous time step. The Laplacian terms are calculated using a standard finite-volume approach. Since the velocity is known at the faces, the gradients can be computed for the cells by a simple application of the Gauss’ Divergence Theorem on the grid cells. Once the gradients are known for the cells, the Gauss’ Divergence theorem is applied again on the dual-grid to obtain a velocity Laplacian at the faces. The face velocity Laplacian at the nodes is computed with a simple inverse distance weighting of the neighboring faces value. Once the Laplacian of the velocity field is known at faces and nodes, the Laplacian term is spatially interpolated using generalized barycentric coordinates to obtain . The location X is the same as in the acceleration term.

Bottom Friction. Bottom friction term is computed semi-implicitly in terms of the bottom friction coefficient, as

5)	$\begin{array}{l}\displaystyle \displaystyle \left ( g S_{f,j} \right )_k = c_{f,k,j}^{n+\theta} V_{k,j}^n+1, \: \: j \in \{ch,lob,rob\}\end{array}$

Bottom friction is computed semi-implicitly where a bottom friction coefficient, , is computed based on -averaged hydraulic radius and velocities. The bottom friction coefficient is therefore updated during the iteration process. At each of those iterations, a new bottom friction term –c_fVⁿ⁺¹ is computed similarly to other implicit terms. The velocity V used in the bottom friction formula is completely implicit for stability purposes.

Solution Procedure. Multiplying the momentum equations for the channel and floodplains by the local area and then summing leads to:

6)

$\begin{array}{l}\displaystyle \displaystyle \sum_{j \in \{ch,lob,rob\} } \frac{ A_{k,j} } { A_k } \left [ \frac{V_{k,j}^{n+1} - V_{X(j)}^n }{\Delta t} +g \left ( \frac{ z_{s,i}^{n+\theta} - z_{s,i-1}^{n+\theta} } { \Delta x_{k,j} } +S_{h,j} \right ) + c_{f,k,j}^{n+\theta} V_{k,j}^{n+1} - \nu_{t,k,j}^n \left ( \frac{\partial^2 V_j } {\partial x^2 } \right )_X^n - \frac{\tau_{s,j}}{\rho h_j} - M_j \right ]\end{array}$

Each term in the above equation is expanded and approximated as

$\begin{array}{l}\sum_{j \in \{ch,lob,rob\} } \frac{A_{k,j}}{A_k} V_{k,j} = V_k\end{array}$
$\begin{array}{l}\sum_{j \in \{ch,lob,rob\} } \frac{A_{k,j}}{A_k} \left ( \frac {V_{k,j} - V_{X(j)} } {\Delta t } \right ) = \frac{V_k - \bar{ V}_X} { \Delta t}\end{array}$
$\begin{array}{l}\bar{ V}_X \approx V_X\end{array}$
$\begin{array}{l}X\end{array}$ = backtracking location corresponding to the velocity $\begin{array}{l}V\end{array}$
$\begin{array}{l}\sum_{j \in \{ch,lob,rob\} } \frac{ z_{s,i}^{n+\theta} - z_{s,i-1}^{n+\theta} } { \Delta x_{k,j} } = \frac{ z_{s,i}^{n+\theta} - z_{s,i-1}^{n+\theta} } { \Delta x_{k,e} }\end{array}$
$\begin{array}{l}X(j)\end{array}$ = backtracking location corresponding to the velocity $\begin{array}{l}V_j\end{array}$
$\begin{array}{l}\sum_{j \in \{ch,lob,rob\} } c_{f,k,j} \approx c_{f,k} = g \frac{A_k^2 } {K_k^2 } |V_k|\end{array}$
$\begin{array}{l}\sum_{j \in \{ch,lob,rob\} } \nu_{t,k,j}^n \left ( \frac{\partial^2 V_j } {\partial x^2 } \right )_X^n \approx \nu_{t,k}^n \left ( \frac{\partial^2 V } {\partial x^2 } \right )_X^n\end{array}$
$\begin{array}{l}\sum_{j \in \{ch,lob,rob\} } \frac{ \tau_{s,j} } {\rho h_j} \approx \frac{\tau_s }{\rho h}\end{array}$

The final semi-discrete 1D momentum equation is given by

7)

$\begin{array}{l}\displaystyle \displaystyle \frac{V_{k}^{n+1} - V_{X}^n }{\Delta t} +g \left ( \frac{ z_{s,i}^{n+\theta} - z_{s,i-1}^{n+\theta} } { \Delta x_{k,e} } +S_{h_k} \right ) + c_{f,k}^{n+\theta} V_{k}^{n+1} - \nu_{t,k}^n \left ( \frac{\partial^2 V } {\partial x^2 } \right )_X^n - \frac{\tau_{s,k}}{\rho h_k}\end{array}$

The approximations made in the backtracking velocity and turbulent mixing term are done for efficiency. The approach avoids computing separate backtracking and interpolations of the velocities and velocity Laplacians for the channel and left and right floodplains.

The momentum equation above can be rearranged to obtain an expression for the velocity at as

8)	$\begin{array}{l}\displaystyle \displaystyle V_{k}^{n+1} = \frac{\Delta t \theta g }{1+\Delta t c_{f,k}^{n+\theta}} \frac{ z_{s,i}^{n+1} - z_{s,i-1}^{n+1} } { \Delta x_{e} } + F_{k}\end{array}$

where

$\begin{array}{l}F_{k} = \frac {B_k^n} { 1 + \Delta t c_{f,k}^{n+\theta} }\end{array}$
$\begin{array}{l}B_k^n = V_X^n - (1-\theta) \Delta t g \frac{z_{s,i}^n - z_{s,i-1}^n}{\Delta x_e} - \Delta t g S_{h,k} + \Delta t \nu_{t,k}^n \left ( \frac{\partial^2 V } {\partial x^2 } \right )_X^n + \Delta t\frac{\tau_{s,k}}{\rho h_k^n}\end{array}$

To obtain a discrete implicit equation for the water volume, the above equation is inserted into the discrete continuity equation to obtain

9)	$\begin{array}{l}\displaystyle \displaystyle \Omega_i^{n+1} + a_{i,i-1} z_{s,i-1}^{n+1} + a_{i,i} z_{s,i}^{n+1} + a_{i,i+1} z_{s,i+1}^{n+1}= b_i\end{array}$

where

$\begin{array}{l}a_{i,i-1} = - \frac{ \Delta t^2 \theta^2 g A_k^{n+\theta} } {(1 + \Delta t c_{f,k}^{n+\theta} ) \Delta x_{e,k}}\end{array}$
$\begin{array}{l}a_{i,i+1} = - \frac{ \Delta t^2 \theta^2 g A_{k+1}^{n+\theta} } {(1 + \Delta t c_{f,k+1}^{n+\theta} ) \Delta x_{e,k+1}}\end{array}$
$\begin{array}{l}a_{i,i} = -a_{i,i-1} - a_{i,i+1}\end{array}$
$\begin{array}{l}b_{i} = \Omega_i^n+\Delta t Q_i + \Delta t \theta \left ( F_k A_k^{n+\theta} - F_{k+1} A_{k+1}^{n+\theta} \right ) + \Delta t (1 - \theta) \left ( V_k^n A_k^{n+\theta} - V_{k+1}^n A_{k+1}^{n+\theta} \right )\end{array}$

There is an equation of this form for every cell in the domain. Before proceeding, the system of equation for all cells is written in a more compact vector notation.

10)	$\begin{array}{l}\displaystyle \displaystyle \boldsymbol{ \Omega ( \boldsymbol{Z} ) } + \boldsymbol{\Psi} \boldsymbol{Z} = \boldsymbol{b}\end{array}$

Where $\begin{array}{l}\boldsymbol{ \Omega }\end{array}$ is the vector of all cell volumes, $\begin{array}{l}\boldsymbol{ Z }\end{array}$ is the vector of all cell water surface elevations, $\begin{array}{l}\boldsymbol{ \Psi }\end{array}$ is the coefficient matrix of the system and $\begin{array}{l}\boldsymbol{ b }\end{array}$ is the right-hand-side vector.

If the coefficients are lagged, the system of equations is mildly non-linear due to the bathymetric relationship for $\begin{array}{l}\boldsymbol{ \Omega }\end{array}$ as a function of $\begin{array}{l}\boldsymbol{ Z }\end{array}$ . The Jacobian (derivative) of Ω with respect to H is given by another bathymetric relationship : the diagonal matrix of cell wet surface areas. If this information is known, a Newton-like technique can be applied to solve the system of equations, producing the iterative formula,

11)	$\begin{array}{l}\displaystyle \displaystyle \boldsymbol{Z}^{m+1} = \boldsymbol{Z}^{m} - \left [ \boldsymbol{\Psi} + \boldsymbol{A}^m \right ]^{-1} \left ( \boldsymbol{ \Omega }^m + \boldsymbol{Z}^m - \boldsymbol{b} \right )\end{array}$

where $\begin{array}{l}m\end{array}$ denotes the iteration index (not to be confused with the time-step).

Robustness and Stability. When there are no fluxes into a cell the water depth is zero, so the water surface elevation is identical to the previous step. In particular, dry cells remain dry until water flows into them.

In the momentum equation, as water depth decreases to zero, all forces tend to zero. However, the bottom friction is the dominant force, so velocities also go to zero in the limit. As a consequence, it is consistent to assume that dry cells have a flow velocity of zero. The momentum equation for dry faces becomes ∂V/∂t=0 and dry faces continue to have zero velocity until water flows into them.

As seen in previous sections, both mass and momentum equations are non-linear. Similarly to the DSW solver, an iterative process must be applied. This idea is presented in steps 6 through 9 below.

The Lagrangian treatment of the of acceleration and mixing terms has the advantage of avoiding stability criteria based on the Courant number and mixing.

Discrete Boundary Conditions. Flow boundary conditions are also discretized:

Water surface elevation: The water surface elevation boundary condition is directly implemented as $\begin{array}{l}z_s^{n+1} = z_{s,b}\end{array}$ . The internal cell is then discretized as described above and the terms containing the boundary water surfaces are placed on the right-hand-side of the system of equations.
Normal Depth: The energy grade slope is specified and utilized to compute a flow at each computation face as $\begin{array}{l}Q^{n+1} = K^n S_{f,b}^{1/2}\end{array}$ . Boundary face flows are included in the internal cells as a source term on the right-hand-side of the system of equations.
Flow: The flow boundary condition is specified at each computational face based on the local conveyance. Boundary face flows are included in the internal cells as a source term on the right-hand-side of the system of equations.

Solution Algorithm. The complete solution algorithm is given here:

The geometry, local orthogonality and subgrid bathymetry data is obtained or pre-computed.
Solution starts with $\begin{array}{l}z_s^0\end{array}$ and $\begin{array}{l}V_k^0\end{array}$ as the provided initial condition at time step n=0.
Boundary conditions are provided for the next time step n+1.
Initial guess $\begin{array}{l}z_s^{n+1} = z_s^{n}\end{array}$ and $\begin{array}{l}V_k^{n+1} = V_k^{n}\end{array}$ .
Compute explicit terms that remain constant through the time step compute, such as the mixing term and wind forcing.
Compute θ–averaged water surface elevation and subgrid bathymetry quantities that depend on the computed elevation (face areas, fluid surface areas, hydraulic radii, Manning’s n, etc.).
A system of equations is assembled for the continuity at the cells as described in the previous sections.
The system of equations is solved iteratively using the Newton-like algorithm with the given boundary conditions to obtain a candidate solution $\begin{array}{l}z_s^{n+1}\end{array}$ .
Velocities $\begin{array}{l}V_k^{n+1}\end{array}$ are computed based on the discrete momentum equation.
If the residual (or alternatively, the correction) is larger than a given tolerance (and the maximum number of iterations has not been reached), go back to step 6; otherwise, continue with the next step.
The computed solution is accepted.
Increment n. If there are more time steps go back to step 3, otherwise end.

The loop provided by steps 6 through 10 has the purpose of updating the coefficients of the system of equations, so that the solution of the nonlinear system (rather than its linearization) is obtained at every time step. As expected, a fully nonlinear solution has very desirable properties such as wetting several cells and propagating waves though several cells in a single time step.