The shallow water equations

In this chapter we study a model for shallow water waves in one dimension:
\begin{align} h_t + (hu)_x & = 0, \label{SW_mass} \\ (hu)_t + \left(hu^2 + \frac{1}{2}gh^2\right)_x & = 0. \label{SW_mom} \end{align}
Here $h$ is the depth, $u$ is the velocity, and $g$ is a constant representing the force of gravity. These equations are "depth averaged" and neglect vertical velocity and any vertical variations in the horizontal velocity. Viscosity and compressibility are neglected, and the pressure is assumed to be hydrostatic. Nevertheless, this is a surprisingly effective model for many applications, particularly when the wavelength is long compared to the fluid depth.

Previously we have looked at a linear hyperbolic system (acoustics) and scalar hyperbolic equations (advection, Burgers' equation, traffic flow). The shallow water system is our first example of a nonlinear hyperbolic system; solutions of the Riemann problem for this system consist of two waves (since it is a system of two equations), each of which may be a shock or rarefaction (since it is nonlinear).

Hyperbolic structure

We can write (\ref{SW_mass})-(\ref{SW_mom}) in the canonical form $q_t + f(q)_x = 0$ if we define
\begin{align} q & = \begin{pmatrix} h \\ hu \end{pmatrix}, & f & = \begin{pmatrix} hu \\ hu^2 + \frac{1}{2}gh^2 \end{pmatrix}. \end{align}
In terms of the conserved quantities, the flux is
\begin{align} f(q) & = \begin{pmatrix} q_2 \\ q_2^2/q_1 + \frac{1}{2}g q_1^2 \end{pmatrix}. \end{align}
Thus the flux Jacobian is
\begin{align} f'(q) & = \begin{pmatrix} 0 & 1 \\ -(q_2/q_1)^2 + g q_1 & 2 q_2/q_1 \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -u^2 + g h & 2 u \end{pmatrix}. \end{align}
Its eigenvalues are
\begin{align} \label{SW:char-speeds} \lambda_1 & = u - \sqrt{gh}, & \lambda_2 & = u + \sqrt{gh}, \end{align}
with corresponding eigenvectors
\begin{align} \label{SW:fjac-evecs} r_1 & = \begin{bmatrix} 1 \\ u-\sqrt{gh} \end{bmatrix} & r_2 & = \begin{bmatrix} 1 \\ u+\sqrt{gh} \end{bmatrix} \end{align}
Notice that -- unlike acoustics, but similar to the LWR traffic model -- the eigenvalues and eigenvectors depend on $q$. Because of this, the waves appearing in the Riemann problem move at different speeds and may be either jump discontinuities (shocks) or smoothly-varying rarefaction waves.

If you wish to examine the Python code for this chapter, see

• exact_solvers/shallow_water.py ... on github,
• exact_solvers/shallow_demos.py ... on github.

For convenience we will work in units such that the gravitational constant $g=1$.

The Riemann problem

Consider the Riemann problem with left and right states

\begin{align*} q_\ell & = \begin{bmatrix} h_\ell \\ h_r u_\ell \end{bmatrix}, & q_r & = \begin{bmatrix} h_r \\ h_r u_r \end{bmatrix}. \end{align*}

Typically the Riemann solution consists of two waves, one related to each of the eigenvectors in (\ref{SW:fjac-evecs}). We refer to the wave related to $r_1$ (the first characteristic family) as a 1-wave, and the wave related to $r_2$ (the second characteristic family) as a 2-wave. Each wave may be a shock or rarefaction. There is a constant intermediate state $q_m$ between them. In the figure below, we illustrate a typical situation in which the 1-wave happens to be a rarefaction and the 2-wave is a shock:

In the figure we have one wave going in each direction, but since the wave speeds depend on $q$ and can have either sign, it is possible to have both waves going left, or both going right. In these cases, the flow is said to be supercritical.

To solve the Riemann problem, we must find $q_m$. To do so we must find a state that can be connected to $q_\ell$ by a 1-shock or 1-rarefaction and to $q_r$ by a 2-shock or 2-rarefaction. We must also ensure that the resulting waves satisfy the entropy condition.

Shock waves

Let us now consider a shock wave connecting a known state $q_*=(h_*, h_* u_*)$ to an unknown state $q=(h,hu)$. In the context of the Riemann problem, $q_*$ will be the known left or right state and $q$ will be the unknown middle state. The Rankine--Hugoniot jump conditions $s(q_* - q) = f(q_*) - f(q)$ for a shock wave moving at speed $s$ read as

\begin{align} s (h_* - h) & = h_* u_* - h u, \\ s (h_* u_* - h u) & = h_* u_*^2 - h u^2 + \frac{g}{2}(h_*^2 - h^2). \end{align}

It is convenient to define $\alpha = h - h_*$. Then we can combine the two equations above to show that the momenta are related by \begin{align} \label{SW:hugoniot-locus} h u & = h_* u_* + \alpha \left[u_* \pm \sqrt{gh_* \left(1+\frac{\alpha}{h_*}\right)\left(1+\frac{\alpha}{2h_*}\right)}\right] \end{align} for a shock moving with speed \begin{align} \label{SW:shock-speed} s & = \frac{h_* u_* - h u}{h_* - h} = u_* \mp \sqrt{gh_* \left(1+\frac{\alpha}{h_*}\right)\left(1+\frac{\alpha}{2h_*}\right)} = u_* \mp \sqrt{gh \frac{h_* + h}{2h_*} }. \end{align}

Choosing the plus sign in (\ref{SW:hugoniot-locus}) (which yields the minus sign in (\ref{SW:shock-speed})) gives a 1-shock. Choosing the minus sign in (\ref{SW:hugoniot-locus}) (which yields the plus sign in (\ref{SW:shock-speed})) gives a 2-shock. Notice from the last expression in (\ref{SW:shock-speed}) that as $h \to h^*$ the shock speed approaches the corresponding characteristic speed (\ref{SW:char-speeds}), while the ratio of the jumps in momentum and depth approaches that suggested by the eigenvectors (\ref{SW:fjac-evecs}).

We can now consider the set of all possible states $(h, h u)$ satisfying (\ref{SW:hugoniot-locus}). This curve in the $h - hu$ plane is known as the Hugoniot locus; there is a family of Hugoniot loci for 1-shocks and another for 2-shocks. Below we plot some members of the 1-loci of curves. In the interactive notebook you can select the 2-loci instead.

The plot above shows the Hugoniot loci in the $h$-$hu$ plane, the natural phase plane in terms of the two conserved quantities. Of course, they all approach the origin as $h \rightarrow 0$. Alternatively, we can plot these same curves in the $h$-$u$ plane:

Note that in this plane the curves in each family are simply vertical translates of one another, and all curves asymptote to $\pm \infty$ as $h \rightarrow 0$. This means that it is impossible to have a shock with an adjacent dry ($h=0$) state.

All-shock Riemann solutions

If we knew a priori that both waves in the Riemann solution were shock waves, we could find $q_m$ by finding the intersection of the 1-Hugoniot locus passing through $q_\ell$ with the 2-Hugoniot locus passing through $q_r$. The Riemann solution would then consist simply of two discontinuities propagating at the speeds given by (\ref{SW:shock-speed}).

Here is an example for which this is the correct solution: The initial depth is the same everywhere, $h_\ell=h_r$, but the velocity $u_\ell>0$ while $u_r<0$ so the flow is toward $x=0$ from both sides. Note that in this plot the dark and light blue stripes show a passive tracer advected with the flow, which is shown only to help visualize the fluid velocity in these plots.

Dam-break problem: all-shock solution

Here's another example, which you might think of as modeling a dam that suddenly breaks. The water is initially at rest on both sides, but much deeper on the left. We can again construct a 2-shock Riemann solution, but this is not physically correct in this case!

Notice that the 1-characteristics (which are plotted as thin lines) don't impinge on the 1-shock; instead, characteristics are diverging away from it. This shock does not satisfy the entropy condition and should be replaced with a rarefaction. The corresponding part of the Hugoniot locus is plotted with a dashed line to remind us that it is unphysical.

Now check the box to plot the 2-characteristics. Notice that these do impinge on the 2-shock; this is an entropy-satisfying shock, and the Hugoniot locus is plotted as a solid line.

The entropy condition

For the 2-shock Riemann solution to be physically correct, each of the waves must satisfy the Lax entropy condition, with the characteristic speed to the left of the wave greater than that to the right. Notice that the shock speed (\ref{SW:shock-speed}) can be written as

\begin{align} s & = u_* \pm \sqrt{gh \frac{h_* + h}{2h_*} }. \end{align}

Thus the 1-shock entropy condition reads as

\begin{align*} u_\ell - \sqrt{gh_\ell} > u_\ell - \sqrt{gh_\ell \frac{h_\ell+h_m}{2h_\ell}} = u_m - \sqrt{gh_m \frac{h_m+h_\ell}{2h_m}} > u_m - \sqrt{gh_m}. \end{align*}

The second expression for the shock speed is obtained by noticing that the shock speed is invariant under transposition of the left and right states. The first and second inequalities both simplify to

\begin{align} \label{SW:1-entropy} h_m > h_\ell. \end{align}

Similarly, it can be shown that the entropy condition for a 2-shock connecting $q_m$ and $q_r$ is satisfied if and only if

\begin{align} \label{SW:2-entropy} h_m > h_r. \end{align}

If (\ref{SW:1-entropy}) or (\ref{SW:2-entropy}) is violated, then the corresponding wave must be a rarefaction rather than a shock.

Go back and look at the phase plane plots above. You'll notice that the 1-Hugoniot locus passing through $q_\ell$ is plotted with a solid line for $h>h_\ell$ and with a dashed line for $h<h_\ell$ (and the same is done for the 2-Hugoniot locus passing through $q_r$). This is to remind us that states connected to $q_\ell$ by the dashed part are unphysical.

Rarefaction waves

Integral curves

In the LWR traffic flow model, we saw that a rarefaction wave consisted of a smoothly varying density. For the shallow water system, both the depth $h$ and the momentum $hu$ will vary smoothly through a rarefaction wave. A given rarefaction wave is associated with just one characteristic family, and the variation within the wave is at all points proportional to the corresponding eigenvector $r_p$. In other words, all values $\tilde{q}=(h,hu)$ in the rarefaction lie on a curve defined by

\begin{align} \label{SW:gen_ode} \tilde{q}'(\xi) & = r_p(\tilde{q}(\xi)). \end{align}

Here $\tilde{q}(\xi)$ is a parameterization of the solution. For the shallow water system, (\ref{SW:gen_ode}) reads as

\begin{align} h'(\xi) = q_1'(\xi) & = 1 \label{SW:dh}, \\ (hu)'(\xi) = q_2'(\xi) & = u \pm \sqrt{gh} = \tilde{q}_2/\tilde{q}_1 - \sqrt{g\tilde{q}_1}, \label{SW:dhu} \end{align}

where we take the minus sign for 1-waves and the plus sign for 2-waves. We can think of (\ref{SW:dh})-(\ref{SW:dhu}) as an initial value ODE; fixing the value of $\tilde{q}$ at one point in the rarefaction wave determines the whole solution of (\ref{SW:dh})-(\ref{SW:dhu}). We refer to each of these solutions as an integral curve.

For the shallow water system, these equations can be integrated exactly. If we fix one point $(h_*, h_*u_*)$ on the curve, the whole integral curve is defined by

\begin{align} hu & = hu_* \pm 2h(\sqrt{gh_*} - \sqrt{gh}), \end{align}

where now the plus sign is for 1-waves and the minus sign for 2-waves. This can be rewritten as

\begin{align} \label{SW:riemann-invariant} u \pm 2 \sqrt{gh} = u_* \pm 2 \sqrt{gh_*}. \end{align}

From (\ref{SW:riemann-invariant}) we see that the value $u + 2 \sqrt{gh}$ is constant along any 1-integral curve, while the value $u-2\sqrt{gh}$ is constant along any 2-integral curve. We refer to these quantities as Riemann invariants for the shallow water system:

\begin{align} w_1(q) & = u + 2 \sqrt{gh}, \\ w_2(q) & = u - 2 \sqrt{gh}. \end{align}

In other words, the trajectories plotted above are just the level curves of $w_1$ and $w_2$. This allows us to easily plot the integral curves as a contour plot of $w_1$ or $w_2$. In the notebook you can select which family of curves to display.

We can also plot the integral curves in the $h$--$u$ plane:

Note that in this plane the integral curves of each family are simply vertical translates of one another due to the form of the functions $w_1$ and $w_2$. Unlike the Hugoniot loci, the integral curves do not asymptote to $\pm\infty$ as $h \rightarrow 0$ and instead each approaches a finite value.

Comparison of integral curves and Hugoniot loci

You may have noticed that the integral curves look very similar to the Hugoniot loci, especially when plotted in the $h-hu$ plane. Below we emphasize the differences by plotting the curve from each of these families that passes through a specified point. The difference is less noticeable in the $h-hu$ plane since all of the curves approach the origin.

In the notebook you can select which family to display.

Near the point of intersection, the curves are very close; indeed, they must be tangent at this point since their direction is parallel to the corresponding eigenvector there (and in fact they also have the same curvature). Far from this point they diverge; for small depths they must diverge greatly, since the Hugoniot locus never reaches $h=0$ at any finite velocity.

The structure of centered rarefaction waves

Rarefactions appearing in the Riemann solution are also similarity solutions; that is, they are constant along any ray $\xi=x/t$. For a $p$-rarefaction connecting states $q_\ell$ and $q_r$, the rarefaction thus has the form

\begin{align} q(x,t) & = \begin{cases} q_\ell & \text{if } x/t \le \lambda_p(q_\ell), \\ \tilde{q}(x/t) & \text{if } \lambda_p(q_\ell) \le \lambda_p(q_r), \\ q_r & \text{if } x/t \ge \lambda_p(q_r). \end{cases} \end{align}

It remains to be determined how $\tilde{q}$ varies as a function of $x/t$. This can be determined by noting that each solution value $\tilde{q}$ propagates at the characteristic speed $\lambda_p(\tilde{q})$. Letting $\xi=x/t$, this gives us the condition

\begin{align} \label{SW:char-sim} \lambda_p(\tilde{q}(\xi)) = \xi. \end{align}

By combining (\ref{SW:char-sim}) with the Riemann invariant condition, we can find $\tilde{h}$ and $\tilde{u}$ within the rarefaction. For concreteness, let us consider a 1-rarefaction connecting $q_\ell$ and $q_m$. Then (\ref{SW:char-sim}) reads $\tilde{u} - \sqrt{g\tilde{h}} = \xi$, which simplifies to

\begin{align} \label{SW:char-sim1} \tilde{h} = \frac{(\tilde{u}-\xi)^2}{g}. \end{align}

Meanwhile, the fact that $w_1(\tilde{q})$ is constant means that

\begin{align} \label{SW:w1-const} \tilde{u} + 2 \sqrt{g\tilde{h}} = u_\ell + 2 \sqrt{gh_\ell}. \end{align}

Combining (\ref{SW:char-sim1}) with (\ref{SW:w1-const}), we find

\begin{align} \label{SW:1-raref-struct-h} \tilde{h} & = \frac{(u_\ell + 2\sqrt{gh_\ell} - \xi)^2}{9g} = \frac{(w_1(q_\ell) - \xi)^2}{9g}, \\ \tilde{u} & = \frac{1}{3}w_1(q_\ell) + \frac{2}{3}\xi. \label{SW:1-raref-struct-u} \end{align}

We see that $h$ varies quadratically through the rarefaction, while $u$ varies linearly.

Connecting states

Suppose we know that the middle state $q_m$ is connected to the left and right states by rarefactions. Then we can find $q_m$ by finding the intersection of the corresponding integral curves. Here is an example for which this is the physically correct solution:

\begin{align} q_\ell & = \begin{bmatrix} 1 \\ -1 \end{bmatrix}, & q_r & = \begin{bmatrix} 1 \\ 1 \end{bmatrix}. \end{align}

In this case the water is flowing away from $x=0$ on both sides.

Notice that the segment of each integral curve that connects to states with a smaller depth is plotted as a solid line, while the segment connecting to states with greater depth is plotted with a dashed line. This again is to remind us that states connected by a rarefaction through the solid part are physical (entropy-satisfying), while states connected by the dashed part would be unphysical (entropy-violating).

Dam-break problem: all-rarefaction solution

Next let us revisit the problem from the section above where we found that the all-shock solution was incorrect. This time we'll try to find an all-rarefaction solution.

Notice that the 2-characteristics (plotted as thin lines) impinge on the 2-rarefaction; in fact they intersect at the left edge of the rarefaction. This means that the solution we constructed is triple-valued and nonsensical as a solution to this one-dimensional conservation law, and so this portion of the solution is omitted in the plots of depth and momentum. In this case a rarefaction wave is not physical and should be replaced with a shock; the corresponding part of the integral curve is hence shown as a dashed line.

The Riemann solution

To determine $q_m$, we need to know whether each of the two waves is a shock or rarefaction, so that we can use the appopriate Hugoniot locus or integral curve. But to determine the nature of each wave, we need to check (\ref{SW:1-entropy}) or (\ref{SW:2-entropy}), which requires knowledge of $h_m$. To deal with this we define a piecewise function $\phi_\ell(h)$ that gives

• for $h>h_\ell$, the momentum connected to $q_\ell$ by a 1-shock;
• for $h<h_\ell$, the momentum connected to $q_\ell$ by a 1-integral curve.

Thus $\phi$ is partly a 1-Hugoniot locus and partly a 1-integral curve; it gives the value of each precisely in the regime where they are physically correct. We can similarly define a function $\phi_r(h)$ related to the 2-Hugoniot locus and 2-integral curve. The middle state depth $h_m$ is the value for which $\phi_\ell(h) = \phi_r(h)$. This value can be found computationally via a rootfinding algorithm.

Dam-break problem: correct solution

Below we plot (at last) the correct solution to the dam-break problem, involving a 1-rarefaction and a 2-shock. Plot the characteristics and notice how each family behaves correctly with respect to the corresponding wave in the Riemann solution.

In the interactive notebook you can confirm that 1-charactersitics spread out across the rarefaction fan while 2-characteristics converge on the shock. View the particle paths and note that between the two waves the fluid velocity is constant, with the fluid accelerating across both the rarefaction and shock to the same intermediate value $u_m$.

Dry states

The approach just described will yield the correct Riemann solution (provided that the rootfinding algorithm converges) in most situations. However, it will fail if the solution involves dry ($h=0$) states.

Dry middle state

Sometimes it seems that there is no way to connect a given pair of left and right states via entropy-satisfying shocks or rarefactions. In the cell below we plot only the physically relevant parts of the integral curves and Hugoniot loci passing through the left and right states. Both the $h-u$ and $h-hu$ planes are shown. With the interactive widget you can try to adjust those states to find a combination for which none of the curves intersect.

You should find that by making the initial states flow sufficiently fast away from each other, there is no intersection in the $h-u$ plane. In the $h-hu$ plane, the integral curves always intersect at the origin. The reason for this ambiguity is that for zero depth it isn't meaningful to assign a velocity. Thus in the $h-u$ plane we could think of the entire $u$-axis as being part of every integral curve. That means that we can always connect the left and right states via an intermediate state with depth $h=0$ (a dry state).

Since the 1-integral curve through $q_\ell$ reaches $h=0$ at $u = u_\ell + 2 \sqrt{gh_\ell}$ and the 2-integral curve reaches $h=0$ at $u=u_r-2\sqrt{gh_r}$, we see that a dry middle state occurs if and only if $u_\ell + 2 \sqrt{gh_\ell} < u_r-2\sqrt{gh_r}$.

It is clear in this case that the solution must consist of two centered rarefaction waves, connecting the left and right states to the middle dry state. These centered rarefactions still have the structure derived above; to complete the solution we only need to know the range of centered characteristics included in each rarefaction. This is most conveniently determined by considering the $h-u$ plane. The 1-integral curve passing through $q_\ell$ is given by $u = u_\ell + 2(\sqrt{gh_\ell}-\sqrt{gh_r})$; it reaches a dry state $h=0$ at $u=u_\ell + 2\sqrt{gh_\ell}$. Similarly, the 2-integral curve passing through $q_r$ reaches a dry state $h=0$ at $u=u_r - 2\sqrt{gh_r}$. Thus the dry state solution is given by

\begin{align} q(x,t) & = \begin{cases} q_\ell & \text{if } x/t \le u_\ell - \sqrt{gh_\ell}, \\ \tilde{q}_\ell & \text{if } u_\ell - \sqrt{gh_\ell} \le x/t \le u_\ell+2\sqrt{gh_\ell}, \\ 0 & \text{if } u_\ell + 2\sqrt{gh_\ell} \le x/t \le u_r - 2\sqrt{gh_r}, \\ \tilde{q}_r & \text{if } u_r - 2\sqrt{gh_r} \le x/t \le u_r+\sqrt{gh_r}, \\ q_r & \text{if } x/t \ge u_r + \sqrt{gh_r}, \end{cases} \end{align}

where $\tilde{q}_\ell$ is given by (\ref{SW:1-raref-struct-h})-(\ref{SW:1-raref-struct-u}) and $\tilde{q}_r$ is given by a similar expression.

Here's an example of a Riemann problem whose solution involves a dry middle state:

Dry initial states

What if the left or right state in the Riemann problem is dry? In this case, the solution consists only of a rarefaction wave connecting the non-dry state to $h=0$.

Interactive phase plane

The live notebook, or this webpage, provides an interactive view that allows you to see how the solution changes as you move the left and right states around in the phase plane. Be sure to click and drag the dark circles indicating the values of $q_\ell$ and $q_r$ in the top left plot. The time can also be adjusted by dragging the circle in the top right plot.

The Hugoniot loci are plotted in red and the integral curves in blue. The portion of each curve that satisfies the entropy condition is plotted as a solid line, while the entropy-violating parts are plotted as dashed lines. The curves are plotted for both the left and right states. Thus the function $\phi_\ell$ just described consists of the solid curve (blue and red parts) passing through the left state $q_\ell$; the function $\phi_r$ is the solid curve (blue and red parts) passing through $q_r$. The intersection of these curves is labeled $q_m$.

Interactive solution

Finally, the notebook contains code that generates an interactive plot showing the solution to the Riemann problem for arbitrary left and right states, including solutions with dry states.