Numerical Hydraulics

2D-Shallow Water Equations

Basics - Solutions - Applications
Cornel Beffa, Dr. sc. techn.

Contents
1

Basics

1.1

Equations and Variables

1.2

Initial and Boundary Conditions

1.3

Validity

Numerical solution

2.1

Standard methods

2.2

Important Terms (explained for 1D case, but analogous in 2D case)

2.3

Accuracy of the cell-centered Finite-Volume-Method (FVM)

2.4

Stability of the cell-centered FVM

2.5

Temporal discretization

Model application

3.1

Grid generation

3.2

Model tests

3.3

Validation

3.4

Areas of application

10

Control questions

11

4.1

Concerning the shallow water equations

11

4.2

Concerning the numerical solution

11

Schwyz, Dec 2008

Dr. C. Beffa, [email protected]

www.fluvial.ch

Hydraulics II Shallow Water Equations: Basics - Solutions - Applications

1 Basics
1.1 Equations and Variables
h = water depth [L]

Variables

q = x- component of specific flow [L2 T-1]

r = y- component of specific flow [L2 T-1]

h q r

0
t x y

Continuityequation

q q 2 gh 2 h

xx

t x h
2

Momentumequations

r qr h

xy
t x h
y

Constants

qr h

z b bx

xy gh
x

by
r 2 gh 2 h

yy gh b

and zb = bottom elevation [L], bi= bottom shear stress [M L-1 T-2], = density of fluid

derived quanti- [M L-3], ii = turbulent normal stresses, ij = turbulent shear stresses

ties
Bottom-shear
stresses

by v u 2 v 2
bx u u 2 v 2

c 2f
c 2f
with (u,v) = x- and y-components of flow velocity and cf = friction coefficient (e.g.
from Manning-Strickler-Formula)

Turbulent shear Influence of turbulent shear stress is important, if velocity gradients are large
stresses

(e.g. in separations, constrictions)

For many applications (e.g. dam break waves, flood waves) the influence of
turbulent stresses can be neglected.
Overview over various approaches in: Rodi W. 1984. Turbulence models and
their application in hydraulics. A state-of-the-art paper. IAHR, Delft

Table 1 - Overview 2D-Shallow Water Equations

1.2 Initial and Boundary Conditions

For the starting time of the simulation at any computational node the values of the 3 variables (h, q, r)
must be given. If no suitable information is available, one of the following options is often chosen in
practice:
(1) Assumption of fluid at rest, constant hydraulic head z = zb +h = const, q = r = 0.

Hydraulics II Shallow Water Equations: Basics - Solutions - Applications

(2) Assumption of a dry channel (dry start): h = q = r = 0.

(3) Using results from an earlier simulation as initial conditions (hot start)
Depending on the flow regime the variables must also be given at the model boundaries. The type of
variables depends on the flow regime and is determined by the number of characteristics which arrive
at the boundary from outside:
Flow regime

Inflow bound.

Ouflow bound.

Subcritical

supercritical

Table 2 Number of required boundary conditions for the 2D-Shallow water equations

Setting meaningful boundary conditions is done in analogy to the 1D case. At an inflow boundary the
flow rate Q should be given as inflow hydrograph together with the direction of the inflow. For an
outflow boundary the following quantities could be used:
(1) Water table elevation (possibly hydrograph)
(2) Water level-flow relationship (if known)
(3) Slope of water surface or energy slope
Suggestion: A comparison with a physical model can simplify the choice of boundary conditions.
Boundaries which can be realized in a physical model will as a rule also lead to reasonable results in a
numerical model.

1.3 Validity
The shallow water equations are based on the assumption that over the flow depth the pressure distribution is hydrostatic. This is the case for long and shallow waves (i.e. waves with a wave length much
larger than water depth), in which the vertical acceleration of fluid elements during the wave passage
stays small. When short or high waves pass, these conditions are violated. In linear wave theory one
distinguishes between long waves (shallow water) and short waves (deep water), which are different
with respect to the movement of water particles (Table below and lecture notes chapter 7).
Region

Shallow water

Transition

Wave length

L/h > 20

2 < L/h < 20

Wave celerity

c gh

Particle path

linear

gL
2h
tanh
2
L

elliptical

Deep water
L/h < 2
c

gL
2

circular

Table 3 Comparison of different linear waves with L = wave length and h = water depth (from Liggett J. A. 1994. Fluid
Mechanics. McGraw-Hill Int. Editions)

Hydraulics II Shallow Water Equations: Basics - Solutions - Applications

During the passage of high waves, significant vertical accelerations can occur, leading to deviation
from the hydrostatic pressure conditions, an essential precondition of the usage of the shallow water
waves. For situations which do not deviate too much from hydrostatic conditions, the classical shallow water equations can be generalized. The Boussinesq-Equations named after the French researcher Boussinesq (1842-1929) contain additional terms which take into account the curvature of
the water surface. They correct the erroneous observation that high waves compared to shallow ones
are slower and subject to larger diffusion.
Under steady state flow the condition of hydrostatic pressure distribution may be violated at discontinuities of the water surface (hydraulic jump, pressure surge) or the channel bottom (steps). The solution of shallow water equations leads to local differences from observations. In such cases there are
also extensions to the shallow water equations, which as a rule lead to new additional equations.
The 2D approach must fail in applications where 3D effects (e.g. secondary flows) become essential.
This happens for strong curvatures of the channel or where the internal turbulent shear stresses influence the flow field considerably. The exact limit up to which a 2D model is still applicable also depends on the accuracy requirements. There are models, which do not average over depth but rather
use the vertical velocity term resulting from the continuity equation. Such 3D shallow water equations
are used for example in the simulation of lakes (e.g. circulation due to wind stress) or in coastal flows.

2 Numerical solution
2.1 Standard methods
The following methods are applied in solving the 2D-shallow water equations: Finite-DifferenceMethod (FDM), Finite-Element-Method (FEM) and Finite-Volume-Method (FVM). Besides the FEM
the following two methods have found a larger fan community:
Staggered Grid Method (SGM)

Cell centered - FVM

Standard method since the 60s

developed in Fluid Dynamics (from

implicit time discretization

1980)

very compact and efficient (Flux computation in explicit time discretisation

parts without interpolation)
usually based on rectangular grids
stable only for subcritical regimes!

requires smooth grids (otherwise inconsistent)

stable for sub- and supercritical regimes

Hydraulics II Shallow Water Equations: Basics - Solutions - Applications

Kontrollvolumen fr
y-Impulsgleichung
Unbekannte

Dreiecksnetz

Kontrollvolumen fr
Kontinuittsgleichung
Kontrollvolumen fr
x-Impulsgleichung

Zellenmittelpunkt

Table 4 Comparison of two standard methods for the solution of the shallow water equations

The staggered grid method has been used since the 60s, above all in coastal flows. A weak point is the
instability for supercritical flows.
The cell-centered FVM is stable for both sub- and supercritical flow regimes. Therefore this method
is the favorite method applied today for inland flows.

2.2 Important Terms (explained for 1D case, but analogous in 2D case)

Consistency: Consistency means that the discrete equations approach (converge to) the differential
equations for x 0 and t 0.
Stability: A stable difference scheme prevents the unlimited growth of numerical error during calculation.
Convergence: The discrete solution Uin approaches the exact solution U(x,t) of the differential equation at every point xi = ix and every time tn = n t, if x 0 and t 0.
Equivalence theorem by Lax: For a correctly posed initial value problem and a consistent discretization, stability is a necessary and sufficient condition for convergence.

2.3 Accuracy of the cell-centered Finite-Volume-Method (FVM)

The cell-centered FVM is of large practical importance due to its stability in both sub- and supercritical flows. To assess the accuracy of the method we look at the one-dimensional advection equation as
a model equation.

u f

0
t x

(1)

with u = vector of unknowns and f = flux vector. A discrete formulation of the cell-centered FVM
with centrally averaged fluxes reads

ui
1
1
x i f i f i 1 f i 1 f i 0
t
2
2

(2)

Hydraulics II Shallow Water Equations: Basics - Solutions - Applications

The accuracy of the discrete formulation can be determined by local Taylor expansion of the fluxes
around point i:
f i 1 f i x

f i x 2 2 f i x 3 3 f i

......
x
2 x 2
6 x 3

(3)

(for nomenclature see figure below).

x-

x+

i-1

i+1

Figure 1 Configuration of computational nodes and no-

xi

menclature

Insertion of expression (3) into (2) and division by xi leads to the differential equation belonging to
discrete equation (2)

u i x x f i
x 2 x 2 2 f i x 3 x 3 3 f i

......
t
2xi
x
4xi
x 2
12xi
x 3

inkonsistent
O ( x )
O x 2

(4)

Comparison of equations (1) and (4) shows:

(1) For equidistant grids (xi = x- = x+) the accuracy of the cell centered-FVM is of second order (Error term O(x) vanishes).
(2) For smooth grids with x = x O(x2) the accuracy is also of second order.
(3) For non-smooth grids the accuracy is only of first order.
(4) For non-smooth grids the discretization is not consistent, as

x x
1.
2 x i

2.4 Stability of the cell-centered FVM

Equation (2) can also be written as

ui
1
x i f i 1 f i 1 0
t
2

(5)

i.e. the flux at point i disappears. A determination of fluxes from the arithmetic average of neighboring values has the consequence that ui becomes independent of fi,. In the results this shows in a numerical solution which is not smooth but has wiggles (compare figure below).
i-1

i+1
i

i+2

Figure 2: Wiggles in the discrete solution when using centrally averaged fluxes

Hydraulics II Shallow Water Equations: Basics - Solutions - Applications

These can be suppressed by so called upwind-methods which introduce an additional numerical diffusion into the scheme. The flux at cell boundary i+ is then
f i 1 / 2

1
f l f r u r u l
2

(6)

where the indices l and r indicate the values to the left and to the right of the cell boundary considered. To determine the factor there are a number of different suggestions, which differ in computational effort and physical motivation (see lecture script chapter 6).
With the obvious choice l = i and r = i+1 for the indices one obtains a scheme of first order accuracy.
For applications, where numerical diffusion is undesirable (e.g. in the simulation of short waves) the
accuracy of the methods can be increased. To do so, the flux terms are determined from the neighboring node values by extrapolation of variables. In order to guarantee a smooth behavior of the solution,
even if there are discontinuities, the flux terms must be limited by so called limiter functions. At the
location of a discontinuity (front) the accuracy is reduced by the limiter to first order.

2.5 Temporal discretization

The discretization in time can be formulated in different ways:

Explicit methods compute the variables of the next time step exclusively from already known
values of the variable at the previous time step. The category of explicit methods also contains the
so called multiple step methods (e.g. predictor-corrector method). Explicit methods are in general
easily programmed but have the restriction that they must fulfill a stability criterion (compare below).

Implicit methods use in addition the values of hitherto unknown values of the variables at the new
time. The solution in point i therefore also depends from the still unknown solution at the
neighboring points of point i. This requires solving of an equations system (as opposed to a single
equation in the explicit method) and due to non-linear terms eventually an iterative procedure.
Implicit methods involve much more effort in programming. At the same time they allow to
choose larger time steps.

Explicit methods gain in importance due to the progress in computer performance. They have to fulfill
the Courant-Friedrich-Levy (CFL) condition. This CFL-number is a measure for the progress of a
disturbance over a time step t related to the grid distance x.
CFL

t
t
u c
x
x

(9)

Hydraulics II Shallow Water Equations: Basics - Solutions - Applications

with u = flow velocity and c = wave velocity. An explicit method is only stable for small CFLnumbers (usually CFL<1.0). Equation (9) thus leads to an equation for the maximum permissible time
step length at fixed grid size x (provided the largest is chosen). In 2D the condition is analogous.

3 Model application
3.1 Grid generation
In the practical application of 2D models, the grid generation is a decisive step. Often it decides upon
success or failure of a model, both with respect to the quality of the results as well as the time consumption for the user.
Application of the FVM requires smooth grids. That means neighboring cells should have similar
dimensions. When using triangles, the appearance of acute-angled triangles (i.e. triangles with one
angle larger 90) should be avoided. Digital terrain models on the basis of TINs (Triangulated Irregular Networks), which are the basis for the spatial discretization in 2D, often contain such triangles.
That means the triangulated DTM is usually not directly utilizable as discretization grid.
A

free

tool

for

the

construction

of

quality

grids

is

the

program

Triangle

(http://www.cs.cmu.edu/~quake/triangle.html ).

Grid generated by buffering of a point raster

Grid generated by quality refinement with Program Triangle

3.2 Model tests

As shown, the FVM are consistent on smooth grids. That means the discretization error decreases
with increasing resolution of the grid. This property can (should) be used to check the model error. To
do so, computations with grids of different refinement are performed and compared. If such solutions
do not show any differences, we have reached grid convergence and the discretization error has no
influence on the results. In the opposite case, the grid resolution has to be increased and the test has to
be repeated.

Hydraulics II Shallow Water Equations: Basics - Solutions - Applications

Figure 3 Model test: Comparison

of the computed water surface
elevations for different grid refinement (Longitudinal profile of the
Danube near Ulm)

3.3 Validation
Calibration: Adjustment of model parameters (e.g. friction coefficients), until results reproduce measurement data sufficiently well.
Inverse Modeling: Automatic calibration of model parameters by adjustment of parameters to reproduce the observed quantities (e.g. Levenberg-Marquardt-Algorithm)
Validation: Calculation of observed quantities, which have not been used in calibration (and without
further change of parameters), and discussion of the differences (including critical discussion of the
quality of the measured quantities).

Provocative hypothesis (Prof. Jean Cunge, University of Grenoble):

By calibration simulation models lose their prognostic capability!

Consequences for the application of flow models (Good Modeling Practice):
(1) Do not overdo the fitting of parameters especially in the light of uncertain observations
(2) Use observations primarily for validation (Learning effect!)

Hydraulics II Shallow Water Equations: Basics - Solutions - Applications

10

3.4 Areas of application

(1) Risk assessment: Simulation of floods, dam and dike break simulation, hydraulic loads of hydraulic structures along and in the river
(2) Hydraulic engineering: Computation of water surface elevations for flows; Design of protective
structures, optimization of structures, abstractions, in- and outflows.
(3) Ecological questions: Assessment of constructive measures on flow (flow velocities, bottom
shear stresses); Assessment of permeability for fish and properties of habitats, water abstractions,
transport of solutes (pollutants) in the far field and intermediate field
(4) Sediment balance: Bottom erosion, flushing, sedimentation of reservoirs

Hydraulics II Shallow Water Equations: Basics - Solutions - Applications

11

4 Control questions
4.1 Concerning the shallow water equations
(1) Which are the unknowns in the shallow water equations?
(2) Which quantities have to be found by empirical formulae?
(3) Give number and types of boundary conditions for inflow and outflow boundaries.
(4) Which assumptions limit the applicability of the shallow water equations?
(5) What is the difference between shallow and deep water waves?
(6) Give typical fields of application for the shallow water equations.

4.2 Concerning the numerical solution

(1) What do the terms consistency, stability and convergence mean in the context of numerical
solution methods?
(2) Under which condition is the Finite-Volume-Method (FVM) consistent?
(3) What is the difference between explicit and implicit solution methods (advantages and disadvantages)?
(4) How is the CFL-number defined and what is its significance in the application of explicit solution
methods?
(5) How are smooth grids characterized and how can they be generated?
(6) How can the size of the discretization error be estimated and how can it be eliminated?