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 can be assembled following a process that mimics the construction of the DSW scheme as

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

where Δ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 n+1). 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. The treatment of the face areas is semi-implicit. This allows for wetting and drying of multiple cells in a single time step and improves the accuracy of the model. However, it can make the solution more difficult and lead to increased iterations.  

Following the same approach used for the DSW equations, the velocities will be expressed as a linear combination of water surface elevation at neighboring cells and terms will be grouped according to their spatial and time indices. All terms related to the time step n will be moved to the right-hand side.

Momentum Conservation

Since velocities are computed on the grid faces, the momentum equations are not located on a computational cell, but rather on a computational face. 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 semi-Lagrangian approach. The Lagrangian form of the advection terms in the momentum equation is computed as:

2) \displaystyle \frac{D \boldsymbol{V}}{Dt} = \frac{\boldsymbol{V}^{n+1}- \boldsymbol{V} ^n_X}{\Delta t}

where \boldsymbol{V} ^{n+1} is located at the computational face and \boldsymbol{V} ^n_X is evaluated at a location X. The backtracking location X is found by integrating the velocity field back in time starting from the location of the computational face. Location X does not in general correspond to a face, so an interpolation technique must be applied to compute \boldsymbol{V} ^n_X. Interpolation in general will not be linear, because the cells are not required to satisfy any constraint in terms of the number of edges. However, it is required that the interpolation algorithm gives consistent values at the boundaries between cells, independently of which cell is used to compute the interpolation. Due to those conditions an interpolation technique based on generalized barycentric coordinates is implemented.

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) g \nabla z_s \cdot \textbf{n}_k = g \frac{\partial z_s^{n+\theta} }{\partial N}


where z_s^{n+\theta} = (1 - \theta) z_s^n + \theta z_s^{n+1}. As a consequence of the implicit weighting, the barotropic pressure gradient consists of two parts, an implicit term with weight \theta and an explicit term with weight (1 - \theta). The water surface Face-Normal Gradient is computed using the neighboring cell water surface elevations. 

Momentum Diffusion

Momentum diffusion represents the sum of molecular and turbulence mixing, and momentum dispersion. Two approaches are available for representing momentum diffusion in HEC-RAS. The first is referred to as the Non-conservative formulation which is based on the Laplacian of the velocity field. The second computes the divergence of diffusive fluxes and is referred to as the Conservative formulation. In HEC-RAS 5.0.7 and earlier, only the non-conservative formulation was available. In general, it is recommended to utilize the Conservative formulation for new models. However, the non-conservative formulation is included for backwards compatibility purposes.

In both formulations, the current velocities are treated explicitly and therefore these terms are computed at the beginning of a time step using the previous time step velocities and water depths.

The momentum diffusion formulation based on the Laplacian of the current velocity is discretized as

4) \left ( \boldsymbol{\nu}_{t} \nabla ^2 u_N \right )_f \approx \nu^n_{t,f} \left ( \nabla^2 \boldsymbol{V} \right )_X^n \cdot \textbf{n}_f

where \nu_{t,k}^n is the explicit face eddy viscosity at face k, \left ( \nabla^2 \boldsymbol{V} \right )_X^n is the Laplacian at location X, and \textbf{n}_k is the face unit vector. 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 \left ( \nabla^2 \boldsymbol{V} \right )_X^n. The location X is the same as in the acceleration term.

The conservative formulation is discretized as

5) \displaystyle \frac{1}{h} \nabla \cdot ( \boldsymbol{\nu}_{t} h \nabla \boldsymbol{V}) \LARGE \vert \normalsize _f \approx \frac{\alpha ^L_f}{\overline{h}_f A_L} \sum_{\begin{array}{c} k \in \text{K}( \text{L}(f)) \\ l \in \text{N} ( \text{L} (f)) \end{array} } A_k \nu_{t,k} \frac{\textbf{n}_f \cdot (\boldsymbol{V}_l - \boldsymbol{V}_L ) }{\Delta x_{L,j}} + \frac{\alpha ^R_f}{\overline{h}_f A_R} \sum_{\begin{array}{c} k \in \text{K}( \text{R}(f)) \\ l \in \text{N} ( \text{R} (f)) \end{array} } A_k \nu_{t,k} \frac{\textbf{n}_f \cdot (\boldsymbol{V}_l - \boldsymbol{V}_R ) }{\Delta x_{R,j}}

in which

              \overline{h}_f = \alpha _f ^L h_L + \alpha _f ^R h_R

              \alpha _f ^L = \frac{\Delta x_f^L}{\Delta x_f^L+\Delta x_f^R}

              \alpha _f ^R = 1 - \alpha _f ^L

               A_k = face vertical area

               A_L = left cell horizontal area

               A_R = right cell horizontal area

               \textbf{V}_j = cell-average current velocity vector

               \nu_{t,k} = \boldsymbol{\nu}_{t,k} \textbf{n}_k \cdot \textbf{n}_{ij}= face eddy viscosity

               \boldsymbol{\nu}_{t,k} = \begin{bmatrix} \nu_{t,xx} & 0 \\ 0 & \nu_{t,yy} \end{bmatrix} = eddy viscosity tensor

               \textbf{n}_k= face-normal unit vector

              \textbf{n}_{ij}  = unit vector for the direction between cell centroids i and j neighboring face k

               \Delta x_{i,j} = distance between cell centroids j and i neighboring face k

The face-normal gradient is approximated by the linear Two-Point Flux Approximation (TPFA) scheme (Edwards and Rogers 1998). The TPFA scheme is robust and monotone. The scheme reduces to the first to second-order central-difference scheme for K-orthogonal meshes. 

The turbulent eddy viscosity is approximated using the longitudinal and transverse components as

6) \nu_{t,k} = \nu_{t,L} + (\nu_{t,L} - \nu_{t,T}) \frac{ ( \textbf{n}_k \cdot \boldsymbol{V}_k )^2 }{\Vert \boldsymbol{V} \Vert_2 }

where

\nu_{t,L} = longitudinal turbulent eddy viscosity [L2/T]
\nu_{t,T}  = transverse turbulent eddy viscosity [L2/T]
\Vert \: \Vert_2 = Euclidean norm operator

Bottom Friction

The nonlinear bottom friction coefficient c_f is computed using the Manning's roughness coefficient as described in Momentum Conservation. As in the previous section, this term will also depend on other quantities, such as the gravitational acceleration, hydraulic radius, Manning’s n and the velocity.

However, extra care must be taken with the bottom friction due to the fact that the term is used implicitly in the equations. Since a Crank-Nicolson type of scheme is being used, the coefficient c_f is computed from -averaged variables located at time and is therefore a \theta-weighted average of the corresponding values at times n and n+1. The bottom friction coefficient c_fis therefore not computed once per time step, but as many times as iterations are required for convergence, through the iteration process of steps 6-10, as it will be seen in the algorithm description below. At each of those iterations, a new bottom friction term -c_fV^{n+1} is computed similarly to other implicit terms. The velocity V used in the bottom friction formula is completely implicit for stability purposes.

Coriolis Effect

The Coriolis term is typically the smallest magnitude term in the momentum equations, but it is also the easiest to compute. The Coriolis parameter f is a pre-computed constant that does not change between time-steps and does not depend on subgrid bathymetry. According to the generalized Crank-Nicolson formula along a streamline, the Coriolis term reduces to Equation 2-179 in the Cartesian oriented system used for the velocities.

7) \displaystyle f_c \boldsymbol{k} \times \boldsymbol{V} \approx \begin{pmatrix} f_c \left[ (1 - \theta) v^n_X + \theta v^{n+1} \right] \\ f_c \left[ (1- \theta) u^n_X + \theta u^{n+1} \right] \end{pmatrix}

The location  where this quantity is interpolated, is obtained from the acceleration advection. As with other implicit terms in the momentum equation, this vector is computed once per iteration.

Fractional Step Method

The solution of the momentum equation uses a fractional-step technique. The first fractional step contains only acceleration and Coriolis terms. The discretization formulas described above yields a vector equation for the velocity. If the coefficients are lagged, this equation is linear on the velocity terms V^n,  and V^{n+1} the water surface elevation terms z_s^n and z_s^{n+1}. The momentum equation contains some velocity cross-terms arising from the Coriolis force. Grouping velocity terms yields the formula:

8) \displaystyle \begin{pmatrix} 1 & \theta f_c \Delta t \\ \theta f_c \Delta t & 1 \end{pmatrix} \begin{pmatrix} u^* \\ v^* \end{pmatrix} = \begin{pmatrix} u^n_X + (1- \theta ) f_c \Delta tv^n_X \\ v^n_X + (1- \theta ) f_c \Delta tu^n_X \end{pmatrix}

where the right-hand side is a linear formula in terms of the velocities and water surface elevations. An explicit formula for V^* without any cross-terms is obtained by: 

9) \displaystyle \boldsymbol{V} ^* = \begin{pmatrix} u^* \\ v^* \end{pmatrix} = \begin{pmatrix} 1 & \theta f_c \Delta t \\ \theta f_c \Delta t & 1 \end{pmatrix} ^{-1} \begin{pmatrix} u^n_X + (1- \theta ) f_c \Delta t v^n_X \\ v^n_X + (1- \theta ) f_c \Delta t u^n_X \end{pmatrix}

The second fractional step adds the acceleration, pressure gradient, eddy viscosity, and bottom friction terms according to the discretization formulas developed earlier.

Solution Procedure

The semi-discrete form of the second fractional step may be written as

10) \displaystyle \frac{D u_{N,f}}{Dt} = \frac{u_{N,f}^{n+1} - u^*_N}{\Delta t} = -g \frac{\partial z_s^{n + \theta}}{\partial N} + M_D - c_{f,f} u_{N,f}^{n+1} + \frac{\tau _{s,N}}{\rho h_f} - \frac{1}{\rho g} \frac{\partial p_a}{\partial N}

where u^* _N = \boldsymbol{V} ^* \cdot \textbf{n} _f is the backtracking velocity including the Coriolis effect from the first fractional step, and the mixing term is computed at the location  X. The term  represents the conservative or non-conservative momentum diffusion term given by 

M_D = \nu^n_{t,f} \left ( \nabla^2 \boldsymbol{V} \right )_X^n \cdot \textbf{n}_f for the non-conservative formulation
M_D = \left[ \frac{1}{h} \nabla \cdot \left( \nu_{t} h \nabla u_N \right) \right]_f^n for the conservative formulation

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

11) \displaystyle u_{N,f}^{n+1} = - \frac{\Delta t g} {1 + \Delta t c_{f,f} } \frac{\partial z_s^{n+\theta}} {\partial N} + F_f

where

F_f = \frac{ B_f^n} {1 + \Delta t c_{f,f}}
B_f^n = u_N^* + \Delta t M_D + \Delta t \left ( \frac{ \tau_s } {\rho h_f^n} - \frac{1}{\rho g} \frac{\partial p_a}{\partial N} \right )

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

12) \displaystyle \Omega_i^{n+1} + \sum_{j \in \text C (i)} a_{i,j} z_{s,j}^{n+1} = b_i

where

a_{i,j} = - \frac{ \Delta t^2 \theta^2 g A_k } { \left ( 1+ \Delta t c_{f,f}^{n+\theta} \right ) \Delta x_N }, 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 - \Delta t \sum_{k \in \text{K} (i) } s_{i,k} \left [ \theta F_k + (1 - \theta) u_{N,k}^n \right] A_k

The system of equations is solved using the same Newton-like iterations used for the DWE solver.

Robustness and Stability

When there are no fluxes into a cell the mass conservation equation becomes h=0, 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.

In contrast with an explicit Eulerian scheme, the semi-Lagrangian scheme for the computation of acceleration has the advantage of avoiding a CFL condition based on velocity.

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 friction slope, S_{f,b},  is specified and utilized to compute a flow at each computation face as Q = K S_{f,b}^{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.
  • The flow boundary condition can be implemented by directly using Equation 2-156. 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 H0 and as the provided initial condition at time step n=0.
  3. Boundary conditions are provided for the next time step n+1.
  4. Initial guess z_s^{n+1}=z_s^n and u_N^{n+1}=u_N^n.
  5. Compute explicit terms that remain constant through the time step compute, such as the diffusion Laplacian field.
  6. Compute θ–averaged water surface elevation and sub-grid bathymetry quantities that depend on the computed elevation (face areas, fluid surface areas, hydraulic radii, Manning’s n, etc.).
  7. A system of equations is assembled for the continuity at the cells as described in the previous sections.
  8. The system of equations is solved iteratively using the Newton-like algorithm with the given boundary conditions to obtain a candidate solution Hn+1.
  9. Velocities u_N^{n+1} are computed based on the discrete momentum equation.
  10. 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.
  11. The computed solution is accepted.
  12. Increment n. If there are more time steps go back to step 3, otherwise end.

The loop provided by steps 6 through10 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.

Similar to the DSW solver, the implementation of this algorithm takes full advantage of computational vectorization and parallelization. Vectorization is used extensively in terms with an explicit discretization in terms of algebraic operations, such as the coefficients for diffusion, Coriolis and bottom friction terms. Simple algebraic steps of the algorithm like steps 4, 6 and 9 are completely vectorized. Parallelization is implemented in terms with a more constructive description such as algorithms, and conditional statements. Examples of such are the subgrid bathymetry table searches, continuity equation coefficients, semi-Lagrangian advection, barotropic pressure gradient terms and the computation of the Laplacian in the eddy diffusion term. Similarly, algorithmic operations in steps 7 and 8 were also parallelized.