Deriving the Shallow Water Equations
This post shows the derivation of the shallow water equations. These equations describe the flow velocity and height of the water’s free surface when the horizontal extent is much greater than the depth. For example, tidal or oceanic flows can be well modelled by these equations. They are derived directly from the Euler equation for inviscid flows.
Contents
Flow Geometry
We consider the depth of the water at each point, with a spatially varying bathymetry. This geometry is shown here:
We can easily convert the fluid depth (\(h\)) and the bathymetry (\(b\)) into a height above (or below) a reference height (\(\eta\)). For geophysical flows, for example, this would be the geoid. Note that \(b\) is the seafloor depth (below reference), so it has the opposite sign to \(\eta\), the height. Consequently we have \(h = \eta - b\).
Derivation
The shallow water equations are derived from the divergence free condition and the incompressible Euler equation, which describe mass- and momentum balance respectively, for the full 3D fluid:
$$\nabla\cdot\mathbf u = 0 \label{div_u}\tag{1}$$$$\Dfrac{\mathbf u}{t} = -\frac{\nabla P}{\rho} + \mathbf f \label{euler}\tag{2}$$
We want to convert these into equations describing the net horizontal components of the velocity, and the fluid depth.
The net horizontal motion of the fluid is obtained by depth averaging the underlying fluid velocity. If we denote the 3D velocity components as \(\mathbf u = (u, v, w)\), then the net velocity’s x-component is given by
$$\bar u(x, y, t) = \frac{1}{h}\zint{u(x, y, z, t)} \label{uavg}\tag{3}$$And similarly for \(\bar v\), meaning we get the 2D velocity field \(\bar{\mathbf u} = (\bar u, \bar v)\).
Conservation of Mass & the Fluid Height
The derivation begins by integrating the divergence free condition for the fluid, \eqref{div_u}, in the vertical direction.
$$ \zint{\Big[{\color{#e2c576} \pfrac{u}{x}} + {\color{#92d8e6} \pfrac{v}{y}} + {\color{#be98df} \pfrac{w}{z}}\Big]}= 0 $$The z term can be handled directly by the Fundamental Theorem of Calculus. For the x and y terms, we want to pull the spatial derivatives outside the integral, but the integrals’ bounds (\(-b\) and \(\eta\)) depend on position. Therefore, in order to pull the spatial derivatives outside, we must take care to apply the Leibniz integration rule correctly. Applying this gives us the following:
$$\begin{align} 0 = &{\color{#e2c576} -u({\small z=\eta}) \pfrac{\eta}{x} - u({\small z=-b})\pfrac{b}{x} + \pfrac{}{x}\zint{u} } \\ &{\color{#92d8e6} -v({\small z=\eta}) \pfrac{\eta}{y} - v({\small z=-b})\pfrac{b}{y} + \pfrac{}{y}\zint{v} } \\ &{\color{#be98df} +w({\small z=\eta}) - w({\small z=-b}) \vphantom{\zint{}}} \end{align}$$Reordering the terms, and using equation \eqref{uavg}
$$\begin{flalign} && 0 = &\pfrac{(h\bar u)}{x} + \pfrac{(h\bar v)}{y} & \\ && &-\Big({\color{#7bd97e} u({\small z=-b})\pfrac{b}{x} + v({\small z=-b})\pfrac{b}{y} + w({\small z=-b}) }\Big) & \text{Seabed} \\ && &-\Big({\color{#d58080} u({\small z=\eta})\pfrac{\eta}{x} + v({\small z=\eta})\pfrac{\eta}{y} - w({\small z=\eta}) }\Big) & \text{Surface} \\ \end{flalign}$$The next step is to identify the boundary condition terms arising from the integrals.
The boundary conditions at both the seabed and the free surface will involve normal vectors. For an arbitrary (differentiable) function \(f(x)\), a vector normal to this function can be found by taking the gradient of a particular scalar field:
$$\mathbf n = \nabla (z - f(x)) = \left[\begin{matrix}-f^\prime(x) \\ 1\end{matrix}\right]$$Note that the resulting vector \(\mathbf n\) does not depend on z. Neither is \(\mathbf n\) a unit vector in general (expect at extrema of \(f\)), however this does not matter for our derivation.
This construction is illustrated here:
Sea-floor Boundary Condition
The first boundary condition is simple enough; the fluid’s velocity at the seabed must be tangential to the seabed itself. If it were not, the fluid would penetrate through the sea floor or leave a vacuum. Mathematically this condition states that the fluid velocity vector is orthogonal to the seabed’s normal vector, i.e.
$${ \mathbf u({\small z = -b}) \cdot \nabla(z + b(x, y))} = 0$$Expanding the gradient and performing the dot-product gives us exactly the seabed boundary terms above (highlighted green), so we see these terms all disappear.
Note that a no-slip boundary condition would give the same result, but is not as general.
Free Surface Boundary Condition
The second boundary condition is slightly subtler. The height of the free-surface at any given position can change in time. However, relative to this vertical motion, the fluid’s velocity must still be tangential to the free surface*. To see why this makes sense, imagine sitting in a reference frame where the free surface (at position \(x\)) is stationary, then we get effectively the same condition as at the seabed.
The vertical speed of the free surface is \(\pfrac{h}{t}\hat{\mathbf z}.\) Thus, using the same construction for a normal vector as above,
$$\Big(\mathbf u({\small z=\eta}) - \pfrac{h}{t}\hat{\mathbf z}\Big)\cdot\nabla(z - \eta(x, y)) = 0$$Expanding this and rearranging yields
$$-\Big({\color{#d58080} u({\small z=\eta})\pfrac{\eta}{x} + v({\small z=\eta})\pfrac{\eta}{y} - w({\small z=\eta}) }\Big) = \pfrac{h}{t}$$So, finally, we can plug this into the equation above for our free surface boundary condition and get
$$\boxed{ \pfrac{h}{t} + \nabla\cdot(h\bar{\mathbf u}) = 0 \label{dhdt}\tag{4} }$$This equation describes mass conservation in the fluid with a free surface. The resulting 2-dimensional flow is no-longer divergence free, we can accommodate non-zero divergence by either “piling up” fluid at a point, or draining it away.
This is the same equation as the usual mass conservation equation for compressible fluids, with the fluid’s height now taking the place of density!
Equation \eqref{dhdt} is in conservation form. By expanding the spatial derivative we can convert it to Lagrangian form:
$$\Dfrac{h}{t} = -h\nabla \cdot \bar{\mathbf u} \label{DhDt}\tag{5} $$The Long-wavelength Limit & Velocity Components
Now we turn to the momentum conservation equations, i.e. Euler’s equation. To make progress here we have to make a long-wavelength assumption.
The long-wavelength limit is that in which things vary over horizontal length-scales much greater than the fluid’s depth. This means the changes in the vertical direction also have to be slow, and the vertical acceleration is approximately zero. Therefore the \(z\) component of the Euler equation becomes
$$-\frac{1}{\rho}\pfrac{P}{z} - g = \Dfrac{w}{t} \approx 0$$where we have used the fact that the vertical component of the body force is gravity alone. Integrating with respect to \(z\) gives
$$P = - \rho g z + C$$We can fix the integration constant by noticing that at the surface we must have atmospheric pressure, so
$$P_0 = - \rho g \eta + C$$thus
$$P = P_0 - \rho g (z - \eta)$$This is the pressure profile we would get in hydrostatic equilibrium. This shows the long-wavelength limit is equivalent to assuming we are approximately at hydrostatic equilibrium. Now we can plug this into the horizontal components of Euler’s equation. Focusing on just the x-component, this gives
$$\begin{align} \Dfrac{u}{t} &= -\frac{1}{\rho}\pfrac{P}{x} + f_x \\[0.2em] &= -g \pfrac{\eta}{x} + f_x \\[0.2em] \end{align}$$All that remains now is to depth-integrate this equation. However, here we run into a snag; the left-hand side is not linear, and so we can’t just pull the integral through the material derivative.
Thus we make our second assumption: that the horizontal components of the velocity do not vary with depth. By expanding the material derivative on the left-hand side, and identifying \(u\equiv\bar u\) and \(v\equiv\bar v\), we get
$$\Dfrac{\bar u}{t} = -g \pfrac{\eta}{x} + f_x \label{DuDt}\tag{6}$$where now the left-hand side is a 2D material derivative, for the 2D depth averaged velocity.
Equation \eqref{DuDt} is in Lagrangian form; let’s convert it to the conservation form.
To do this it is easier to consider the total x-momentum \(h\bar u\)1, rather than the velocity \(\bar u\). The material derivative obeys the usual product rule, therefore
$$\begin{align} \Dfrac{(h\bar u)}{t} = \pfrac{(h\bar u)}{t} + \bar{\mathbf u}\cdot \nabla (h\bar u) = h\Dfrac{\bar u}{t} + \bar u\Dfrac{h}{t} \\[0.2em] \end{align}$$Now plug in equations \eqref{DhDt} and \eqref{DuDt} for the material derivatives on the right-hand side, and use the product-rule in reverse for the spatial gradients (highlighted in yellow)
$$\begin{align} \pfrac{(h\bar u)}{t} + {\color{#e2c576} \bar{\mathbf u}\cdot \nabla (h\bar u)} &= h\Big(-g\pfrac{\eta}{x} + f_x\Big) + {\color{#e2c576} \bar u (-h\nabla\cdot \bar{\mathbf u})} \\[0.2em] \pfrac{(h\bar u)}{t} + {\color{#e2c576} \nabla\cdot (h\bar u \bar{\mathbf u})} &= -gh\pfrac{\eta}{x} + hf_x \\[0.2em] \end{align}$$The derivation is identical for the y component.
This completes the derivation. Henceforth, we will only ever consider the two dimensional velocity, so I’ll drop the overbar.
Shallow Water Equations in Conservation Form
The full shallow water equations in conservative form read
$$\begin{align} \pfrac{h}{t} + \nabla \cdot (h \mathbf u) &= 0 \\[0.2em] \pfrac{(h\mathbf u)}{t} + \nabla \cdot (h \mathbf u \otimes \mathbf u) &= -gh \nabla \eta + h \mathbf f \end{align}$$The first equation has no source terms (zero right-hand side), which tells us that mass is exactly conserved. The water column depth can (and does) change both in space and time, but throughout the system as a whole, fluid is neither created nor destroyed.
The second equation does have source terms, telling us that momentum is not conserved. This is not surprising, the change in momentum comes exactly from the body forces experienced by the fluid. These being gravity for the vertical direction, and any other forces in the horizontal directions (e.g. Coriolis force.)
It is a simple matter to write these in terms of the displacement height \(\eta\), we just substitute in \(h= \eta + b\), yielding
$$\begin{align} \pfrac{\eta}{t} + \nabla \cdot (\eta \mathbf u) &= {\color{#92d8e6} -\nabla\cdot(b\mathbf u)} \\[0.2em] \pfrac{(\eta\mathbf u)}{t} + \nabla \cdot (\eta\mathbf u \otimes \mathbf u) &= -gh \nabla \eta + h \mathbf f \,{\color{#92d8e6} - \nabla\cdot(b \mathbf u \otimes \mathbf u)} \end{align}$$We gain extra source terms for both equations (highlighted in blue), which represent the effect of the bathymetry on the fluid’s height.
Comments on the Derivation
In order not to break up the flow of the derivation with asides and notate bene, I have put these observations into this separate section.
Tong presents a very terse & understandable description of the derivation, but with the caveat that he treats only the case of uniform bathymetry. Because the seabed has a uniform depth, the bottom boundary condition becomes a condition purely on \(w,\) and the Leibnitz terms arising from moving the depth integral inside the horizontal spatial derivatives become zero.
Free Surface Boundary Condition
The free surface boundary condition seems to be presented differently between sources. For example, in Segur and Tong it is given as
$$\Dfrac{\eta}{t} = w({\small z = \eta})$$Rather than the “no relative normal flow” condition used in Dawson & Mirabito.
This seems to make sense, but it actually relies on identifying the x and y components of the 3D velocity at the surface, with the 2D depth averaged velocity. This is something you can validly do only if you’ve already made the second assumption (uniform horizontal velocity profile).
Only the “no relative normal flow” criterion correctly applies to the 3D velocity, avoiding the need to make the uniform horizontal velocity profile at this early point in the derivation. (In Segur they have not explicitly made this assumption when using the boundary condition. However, they do make it later on, meaning it comes out in the wash.)
One doubt I had about the “no relative normal flow” constraint was; why we do not need to account for the fact the free surface is changing shape? I suppose it’s because the divergence free condition is something that must be obeyed at all instants of time and at every location separately.
Depth Integration of the Euler Equation
At the point we integrate the Euler equation we made the uniform horizontal velocity profile assumption. This lets us identify the horizontal components of the 3D velocity with the 2D depth averaged velocity, and side-step the vertical integration.
Had we not made this assumption, we would have obtained an equation very similar to equation \eqref{DuDt}, but with extra terms arising from the interaction of the non-linear advection term and the integration. Dawson & Mirabito mention this, but they do not give the form of these terms. They refer to them simply as “differential advection terms” (differential here meaning the velocity is different at different depths,) and make the second assumption by assuming these terms are zero.
Analogous to the way the distribution of microscopic velocities in the Boltzmann equation leads to diffusion of momentum (see here), I expect these terms (which are associated with the distribution of horizontal velocities) would manifest as a kind of “viscosity”.
Linearization & Gravity Waves
In their full glory, the shallow water equations contain many non-linear terms. We can linearize them by assuming the velocity \(\mathbf u\) and height \(\eta\) are only small perturbations away from zero. In this way, any cross terms of \(\mathbf u\) and \(\eta\) are even smaller, and we can neglect them. So, simply dropping all non-linear terms in equations \eqref{dhdt} and \eqref{DuDt} gives
$$\begin{align} \pfrac{\eta}{t} &= - \nabla \cdot(b\mathbf u) \label{detadt_linear}\tag{7} \\[0.5em] \pfrac{\mathbf u}{t} &= -g\nabla \eta + \mathbf f \label{dudt_linear}\tag{8} \end{align}$$Note that the water column depth \(h=\eta+b\) may be large because \(b\) can be large. Moreover, remember that \(b\) is a fixed profile, so it does not count as making things non-linear.
It is simple to show that the linear Shallow water equations give rise to gravity waves, as one might expect. For this, let’s assume a constant depth \(b=B\), and no horizontal body force (\(\mathbf f = 0\)). First, take the time derivative of equation \eqref{detadt_linear} then substitute in equation \eqref{dudt_linear}:
$$\begin{flalign} && \partial_{tt}{\eta} &= -\partial_t \nabla\cdot(B\mathbf u) & \text{Differentiate} \\[0.5em] && &= -B\nabla\cdot (\, \partial_t\mathbf u) & \text{Commutativity} \\[0.5em] && &= -B\nabla\cdot (-g\nabla \eta) & \text{Substitute \eqref{dudt_linear}} \\[0.5em] && &= gB\nabla^2 \eta & \text{Distribute} \end{flalign}$$which is the wave equation with \(c = \sqrt{gB}\). This matches the result one obtains from dimensional analysis.
This derivation is similar to how we show that electromagnetic waves arise from Maxwell’s equations.
The fact the wave speed depends on depth is often used as an explanation for why waves break on the beach. The argument goes like this: as a wave climbs the beach the depth, and thus \(c\), decreases. Thus the faster water from behind catches up to the slower water ahead. Eventually there is so much water piling up that it becomes unstable and overtops.
If we don’t assume constant depth and no body force, we merely gain a couple of extra terms, turning our wave equation into an inhomogeneous wave equation.
Alternative Derivation
Here is another derivation of the shallow water equations, which utilizes the divergence theorem. This derivation is based directly on the notion of conservation, and so it directly gives us the conservative form.
Consider a finite sized water column, of arbitrary cross-section. Denote the cross-section of the water column by \(A\), its boundary by \(\partial A\), and the unit normal to the boundary as \(\hat{\mathbf n}.\)
The derivation considers the flow of fluid & momentum through the sides, and into the water column.
Mass Conservation
The inflow of fluid through any side is given by
$$q = \int_{-b}^\eta \mathbf u\cdot(-\hat{\mathbf n})\,dz = -(h\bar{\mathbf u})\cdot\hat{\mathbf n}$$\(\hat{\mathbf n}\) is the outward unit normal, hence the minus sign. We can write the net inflow as
$$Q = - \oint_{\partial A} (h\bar{\mathbf u})\cdot\hat{\mathbf n} \,dl $$The divergence theorem tells us
$$Q = - \int_A \nabla\cdot (h\bar{\mathbf u})\,dA$$Now coming from the other side, because our fluid is incompressible, the change in volume of fluid inside our region \(A\) must manifest as a change in the water’s depth. Therefore
$$Q = \int_A \pfrac{h}{t} \,dA$$We did not specify anything about the size nor shape of region \(A\), therefore the integrands themselves must be equal giving
$$\pfrac{h}{t} = -\nabla\cdot (h\bar{\mathbf u})$$which is the same as our mass conservation equation \eqref{dhdt}.
I actually find this derivation a bit more natural, but I have not seen it presented anywhere. Perhaps it contains an error I’ve not spotted? Maybe using the incompressibility to say flow results in a height change, while physically correct, is jumping the gun?
Momentum Conservation
For the momentum, this conservation view point is not simpler than the derivation presented above, because there identifying \(u\equiv\bar u\) avoids the depth integration. However, I will include the alternative derivation, partly for completeness, and partly because it gives a different perspective (Eulerian vs Lagrangian.)
As above, I’ll show the derivation for the x-component of velocity only. The y-component is derived in exactly the same manner.
The x-momentum flow into the water column through any side is given by
$$\begin{align} j = - \zint{u\mathbf u \cdot \hat{\mathbf n}} = -(h\bar u \bar{\mathbf u} + \Delta\bar{\mathbf u})\cdot\hat{\mathbf n} \end{align}$$where \(\Delta\bar{\mathbf u}\) represents the “differential advection terms” discussed above. From here on we make the uniform horizontal velocity profile assumption and set \(\Delta\bar{\mathbf u} = 0.\)
We then integrate this over the whole surface of the water column to get the total momentum flux into the volume, and apply the divergence theorem
$$\begin{align} J &= - \oint_{\partial A} h\bar u\bar{\mathbf u}\cdot\hat{\mathbf n} \,dl \\[0.2em] &= - \int_{A} \nabla\cdot( h\bar u \bar{\mathbf u}) \,dA \label{J_as_div}\tag{9} \end{align}$$Similarly to above, let’s consider the change of momentum throughout the volume. Here, however, we must be careful because we now have sources of momentum in the volume; namely the forces the fluid feels. The total rate of change of momentum at any point is the sum of the advection (momentum flow) and source (force) terms. Thus the change of x-momentum in our volume due only to advection is
$$J = \int_{A}\zint{\pfrac{u}{t} - \Big(\!-\!\frac{1}{\rho}\pfrac{P}{x} + f_x\Big)}\,dA$$Now we make the long-wavelength / hydrostatic assumption exactly as above, then vertically integrate.
$$\begin{align} J &= \int_{A}\zint{\pfrac{u}{t} - \Big(\!-\!g\pfrac{\eta}{x} + f_x\Big)}\,dA \\[0.2em] &= \int_{A} \pfrac{(h\bar u)}{t} - \Big(\!-\! hg \pfrac{\eta}{x} + hf_x\Big)\,dA \\[0.2em] \end{align}$$We have again used the identification \(u\equiv\bar u\) in this step.
Equating this with equation \eqref{J_as_div} and noting that \(A\) was completely arbitrary gives immediately
$$\pfrac{(h\bar u)}{t} + \nabla\cdot(h\bar u\bar{\mathbf u}) = -hg\pfrac{\eta}{x} + hf_x$$completing the derivation.
References
-
We can ignore density since it cancels everywhere at the end. ↩︎