Discrete Scheme for DWE

An overview of the discretization and solution algorithm of the DWE is described in section. The derivation begins with the discretization of continuity equation and momentum equations. The discrete form of the SWE is then obtained from the continuity and momentum equations. Finally, the solution algorithm is described. 

Continuity Equation

The continuity equation is discretized for the DWE solver as:

1) \displaystyle \frac{ \Omega _i ^{n+1} -\Omega _i ^n} {\Delta t } + \sum_{k \in \text{K}(i)} s_{i,k} u_{N,k}^{n+\theta} A_k^n = Q_i

where \Delta t is the time step, and the velocities have been interpolated in time using the generalized Crank-Nicolson method (which is used to weight the contribution of velocities at time steps  and ). Since the momentum equation is rotation invariant, it will be assumed that s_{i,k} is the sign in the outward direction at face k. Note that the face areas are treated explicitly. This increases the stability of the solver but reduces the accuracy for large time steps and limits the wetting and drying to one cell at a time for each time step.

Momentum Equation

The diffusive wave approximation to the momentum equation can be written in discrete form as

2) \displaystyle u_{N,k}^{n+\theta} = -\alpha_k \frac{z_{s,R}^{n+\theta} -z_{s,L}^{n+\theta}}{\Delta x_{N,k}}

where

\alpha_k = \frac{ R_k^{2/3}} {n | \nabla z_s |^{1/2} }

z_s^{n+\theta} = (1 - \theta) z_s^n + \theta z_s^{n+1}
\Delta x_{N ,k} = face-normal distance between cells L and R 

Diffusion-Wave Equation

The DWE describes the conservation of mass and momentum. A discrete DWE is obtained by substituting the diffusion-wave approximation for the velocity into the above continuity equation leading to

3) \displaystyle \frac{ \Omega _i ^{n+1} -\Omega _i ^n} {\Delta t } - \sum_{k \in \text{K}(i)} s_{i,k} \frac{ \alpha_k} {\Delta x_N } \left ( z_{s,R}^{n+ \theta} - z_{s,L}^{n+ \theta}\right ) = Q_i

Once the DSW equation has been solved, the velocities can be recovered by substituting the water surface elevation back into the Diffusion Wave equation. The above equation may be written in compact form as

4) \displaystyle \Omega _i ^{n+1} + \sum_{j\in\text{C}(i)} a_{i,j} z_{s,j} ^{n+ \theta} = b_i

where

a_{i,j} = - \frac{ \Delta t \theta \alpha_k } {\Delta x_{N,k} }, \forall \: j\in\text{N}(i)

a_{i,i} = - \sum_{j\in\text{N}(i)} a_{i,j}

b_i = \Omega_i^n + \Delta t Q_i - \frac{1-\theta} {\theta} \sum_{j\in\text{C} (i)} a_{i,j} z_{s,j}^n

The system of equation for all cells is written in a more compact vector form as

5) \boldsymbol{\Omega} (\boldsymbol{Z}) + \boldsymbol{\Psi} \boldsymbol{Z} = \textbf{b}

Where \boldsymbol{\Omega} is the vector of all cell volumes, \boldsymbol{Z} is the vector of all cell water surface elevations at time n+1, \boldsymbol{\Psi} is the coefficient matrix of the system and \textbf{b} is the right-hand-side vector.
The system of equations is mildly non-linear due to the bathymetric relationship for \Omega as a function of \boldsymbol{Z} . The Jacobian (derivative) of Ω with respect to \boldsymbol{Z} is given by another bathymetric relationship \boldsymbol { A(z_s)} : 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,

6) \boldsymbol{Z}^{m+1} = \boldsymbol{Z}^m - \left[ \boldsymbol{\Psi} + \boldsymbol{A}^m \right] ^{-1} \left( \boldsymbol{\Omega^m} + \boldsymbol{\Psi} \boldsymbol{Z}^m - \boldsymbol{b} \right)

where m denotes the iteration index (not to be confused with the time-step).

Robustness and Stability

When there are no fluxes the coefficients are a_{i,j}=0 and Q=0 so the water surface elevation is identical to the previous step. In particular, dry cells remain dry until water flows into them. Equation (4) implies that the coefficients a_{i,j} of Equation 2-170 will depend on the value of the water surface elevation. In order to maintain consistency with the generalized Crank Nicolson method, the terms a_{i,j} must be evaluated at time n+ \theta, creating a circular dependence on the solution of the system of equations. This situation is typical of nonlinear systems and is corrected through iteration. This is presented in steps 5–8 below.
The linearized scheme is unconditionally stable for 1/2 \leq \theta \leq 1. When \theta \leq 1/2the scheme is stable if:

7) \displaystyle \frac{\Delta t}{(\Delta x)^2} < \frac{1}{2 - 4 \theta}


When \theta = 1, the scheme obtained is implicit. It corresponds to using backward differences in time and positioning the spatial derivatives at step n+1. When \theta = 1/2, this is the Crank-Nicolson scheme obtained from central differences in time and positioning the spatial derivatives at n+1/2.
The linearized scheme is second order accurate in space. The time accuracy depends on the choice of \theta; for instance, for \theta =1 it is first order accurate and for \theta =1/2 it is second order accurate.

Discrete Boundary Conditions

Flow boundary conditions are also discretized:

  • Water Surface Elevation: The water surface elevation boundary condition is directly implemented as z_s^{n+1} = z_{s,b}. 
  • Normal Depth: The energy grade slope is specified and utilized to compute a flow at each computation face as Q = K S_f^{1/2} . 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 similarly implemented as a condition on the water surface gradient using a finite volume approximation of Manning's equation. A rotation of the local coordinate system is necessary since, in general, the direction normal to the boundary does not coincide with the Cartesian directions.

Solution Algorithm

The complete solution algorithm is given here:

  1. The geometry, local orthogonality and sub-grid bathymetry data is obtained or pre-computed.
  2. Solution starts with z_s^0 as the provided initial condition at time-step n = 0.
  3. Boundary conditions are provided for the next time step n+1 .
  4. Initial guessz_s^{n+1} = z_s^n.
  5. Compute the \theta-averaged water surface elevation z_s^{n+\theta} and subgrid bathymetry quantities that depend on it (face area, horizontal surface area, hydraulic radius, Manning's n, etc.).
  6. The coefficients a_{i,j} are computed and the system of equations is assembled.
  7. The system of equations is solved iteratively using the Newton-like algorithm with the given boundary conditions to obtain a candidate solution z_s^{n+1}.
  8. If the residual (or alternatively, the correction) is larger than a given tolerance (and the maximum number of iterations has not been reached), return to step 5; otherwise continue with step 9.
  9. The computed z_s^{n+1} is accepted and the velocities u_{N}^{n+1} can be calculated using the discrete version of the momentum equation.
  10. Increment n. If there are more time steps go back to step 3, otherwise end.

The loop provided by steps 5 through 8 has the purpose of updating the coefficients a_{i,j} 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 or updating coefficients that are evaluated at time n+ \theta.