The Shallow Water Equations

The shallow water equations are frequently used for modeling both oceanographic and atmospheric fluid flow. Models of such systems lead to the prediction of areas eventually affected by pollution, coastal erosion, and polar ice-cap melting.

Comprehensive modeling of such phenomena using physical descriptions such as the Navier-Stokes equations can often be problematic, due to the scale of the modeling domains as well as through resolving free surfaces. The shallow water equations, of which there are a number of representations, provide an easier description of such phenomena.

This 1D model investigates the settling of a wave over a variable bed as a function of time. The initial wave and the shape of the bed are represented by mathematical relations so that it is easy to change parameters such as the wave amplitude or the bed’s shape.

This example uses the Saint-Venant’s shallow water equations, which are the following:

and

where z is the thickness of the water layer (m), v is the velocity (m/s), g is the gravity constant (m/s2), and ν is the kinematic viscosity (m2/s). The definition of the thickness of the water layer, z, is zs − zf, where zs and zf are the measures in Figure 1 below. For further details, see Ref. 1.

Figure 1: Representative vertical section through the fluid domain showing the bed of a lake and the water surface.

Artificial Stabilization

With time, the flow develops discontinuities known as hydraulic jumps. Use artificial stabilization to replace the jumps by steep fronts that can be resolved on the grid. Small amplitude waves on still water of depth z move with velocity

. The maximal propagation velocity is

for water waves.

Stabilize the solution by adding artificial viscosity chosen to make the cell Reynolds number of order unity. To do so, add the term tune vphase h to the kinematic viscosity ν. Here, tune is an O(1) tuning parameter and h is the local element size. You add the contribution to the divergence term of the conservation law so that it does not affect the shock speeds. The modification is first order in element size.

Model Definition

This application studies a simple example of shallow water in a channel with bottom topography shown in Figure 1. Notice the difference in scale between the x and y directions.

Figure 2: Sea bed profile, zf, used in the model.

Constraints (v = 0) are implemented at both ends, while the physics are described by the equations above. The initial condition is a wave profile, which the following expression defines:

where zf is the analytical expression for the sea bed profile (see Figure 2). The elevation of the water surface is z + zf, while Figure 3 shows z0 + zf.

Figure 3: The initial water surface profile, z0 + zf, and the sea bed profile, zf.

Results and Discussion

The simulation runs for 60 seconds. Figure 4 shows the water surface and slope of the sea bed at six output times toward the beginning of the simulation.

Figure 4: The water level and the slope of the sea bed at six output times. Time spans from t = 0 to t = 15 at steps of 3 seconds.

The simulation clearly shows the influence of the topography of the sea bed on the elevation of the water surface. Another interesting visualization of the results is an animation, which is easy to create using COMSOL Multiphysics.

Notes About the COMSOL Implementation

The modeling procedure is straightforward using the General Form PDE interface with two dependent variables: Z and V. It is easy to define expressions, such as the one that describes the initial wave profile, z0, as a variable in the model.

Reference

1. O. Pironneau, Finite Element Methods for Fluids, John Wiley & Sons, 1989.

COMSOL_Multiphysics/Equation_Based/shallow_water_equations

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the table, enter the following settings:


Z
V

Click

In the tree, select .

Click

Root

In the window, click the root node.

In the root node’s window, locate the section.

From the list, choose .

Because the dependent variables, Z and V, are of different dimensions, it is convenient to switch off unit support and keep track of the units by hand instead.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
nu1	1e-6	1E-6	Kinematic viscosity (m^2/s)
x0	6	6	Sea bed ridge position (m)
a	0.005	0.005	Sea bed ridge height (m)
k1	0.0015	0.0015	Sea bed slope parameter
tune	0.1	0.1	Tuning parameter

Geometry 1

Interval 1 (i1)

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Coordinates
0
10

Click

Definitions

Variables 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Description
Zf	aexp(-(x-x0)^2)+k1x	Sea bed profile
dZfdx	d(Zf,x)	Sea bed profile, x derivative
Zs	Z+Zf	Sea surface level
Z0	0.02-Zf+0.005*exp(-(x-3)^2)	Initial water depth profile
vphase	abs(V)+sqrt(g_const*Z)	Maximal wave propagation velocity
nu	nu1+vphasehtune	Effective kinematic viscosity