Surprisingly Rich in Applications: The Shallow Water Equations and Numerical Solution Methods

By Jörn Behrens

When I started my scientific career in atmospheric modeling, shallow water equations were considered a simple equation set; complex enough to serve as a useful test set, but not more than that [13]. Then I changed my application field to tsunami forecasting and found that shallow water equations have much more potential than my previous experience indicated. Therefore, a more careful view of this seemingly simple (yet interestingly complex) nonlinear equation set appears worthwhile. The eventual motivation for this article was triggered by my participation in a minisymposium entitled “New Numerical Approaches to Ocean Modeling Across Scales,” held at the 2016 SIAM Annual Meeting in Boston, Mass., last summer.

Equations

There are several ways to motivate this system of equations. Currently the system is often derived from Navier-Stokes equations via integration over the fluid depth. The argument is that a shallow fluid stratum (characterized by typical ratios \(L\gg H\) and \(H\gg h\), where \(L\) is typical horizontal length scale, \(H\) is mean water depth, and \(h\) is typical water level elevation over the mean) has only negligible vertical variability and can thus be modeled in a homogeneous two-dimensional layer.

Another approach derives the equations from first principles, considering mass and momentum conservation with Newton’s first and second laws in mind. French engineer and mathematician A.J.C. Barré de Saint-Venant, who was interested in hydraulic engineering and open channel flow, first discovered this approach in a one-dimensional version during the second half of the seventeenth century [11]. As a result, shallow water equations are also referred to as Saint-Venant equations.

This nonlinear coupled system of two basic \(d\)-dimensional equations is given by the continuity and the momentum equations:

\[\frac{\partial h}{\partial t} + \nabla\cdot(\mathbf{v}(H+h)) = 0, \tag1 \]

\[\frac{\partial \mathbf{v}h}{\partial t} + \nabla\cdot(\mathbf{v}\ \mathbf{v}h) + g\nabla h = s(\mathbf{v},h), \tag2\]

defined in a space-time domain \(\tilde{\Omega}\subset\mathbb{R}\times\mathbb{R}^d\) with suitable initial and boundary conditions. We denote the free surface elevation with \(h\), mean fluid layer depth with \(H\), a \(d\)-dimensional velocity vector with \(\mathbf{v}\), the gravitational acceleration with \(g\), and a function of sources with \(s\), which can be among other frictional forces, internal viscosity, or stress terms. \((1)\) and \((2)\) are coupled over the unknown prognostic variables \(h\) and \(\mathbf{v}\). While \(h\) describes the propagation of water elevation by velocity (and constitutes the conservation of mass), \(\mathbf{v}\) establishes the propagation of momentum subject to external and internal forces.

Applications

It thus becomes quite obvious that shallow water equations are well-suited for wave propagation phenomena. In fact, we can show that the set of two first-order equations given in \((1)\) and \((2)\) is the counterpart of the second-order wave equation. So, as previously mentioned, this equation type has shown great benefit in fields such as tsunami modeling [2]. Most—if not all—of the current operational tsunami warning centers worldwide rely on shallow water theory-based modeling. While we know that dispersion and non-hydrostatic pressure affect the wave behavior in certain situations, the uncertainty in determining initial conditions in case of a tsunami event dominates the overall forecasting accuracy. Robustness and relative simplicity of the two-dimensional equation set provides valuable information for assessing the hazard.

Many application fields benefit from shallow water theory. From wave propagation in tsunami forecasting [3] over storm surge simulations [4] to river flooding and dam break problems [9], this type of equation set successfully treats a large number of hydrodynamical issues. Shallow water equations also model landslides and avalanches [5]. Of course, the limitations are quite restrictive in this field of application, and validity of results must be carefully evaluated.

The Advanced Simulation of Complex Earthquake Tsunami Events (ASCETE) project¹ performs high-fidelity earthquake-tsunami simulations to gain better knowledge about the interaction of complex rupture mechanisms with wave propagation and inundation behavior. Figure 1 offers one such example, where a yellow-red colormap shows the non-uniform rupture in terms of uplift. A purple-blue colormap for the ocean surface displacement/wave height shows the corresponding tsunami wave.

Figure 1. Complex rupture-tsunami simulation: earthquake generated non-uniform uplift (red-yellow color map) and sea surface displacement (blue-purple color map).

Solution Methods

One interesting feature of shallow water theory is its non-linearity. This equation set represents one of the simplest, however relevant, settings in which one can study nonlinear effects. While this is interesting, the observation of surprising interactions (e.g., nonlinear steepening of tsunami waves in interaction with tides [1]) poses severe challenges to the numerical method for solving them. Therefore, researchers have presented a vast number of different approaches. Besides traditional finite difference type discretization methods, finite volume methods are well-established and most efficient when combined with adaptive mesh refinement [8, 10]. Finite element methods are also popular because they can accommodate for arbitrary cell shapes and are thus considered well-suited for complex domain shapes and topographical features [6].

Delicate balances often drive the main dynamics of geophysical application fields. Therefore, well-balanced numerical methods are of great importance. One can explain well-balancedness in a simple way: when setting velocity to zero in equation (2)—this represents a fluid in rest—one can see that the gradient term and the source term must cancel. This requirement is fulfilled naturally in continuous form. However, it is a non-trivial requirement for many numerical schemes and crucial for accurate results [12].

Future

While the application fields are wide and results of shallow water theory often very satisfying, severe restrictions naturally exist. In particular, when three-dimensional effects like flow over steep topography, wave dispersion, and non-hydrostatic vertical pressure variations play a role, the simple two-dimensional model breaks down. 3D modeling is computationally too demanding in many cases, which is why one instead considers extensions of the shallow water equations. Recently, we were able to show that a commonly-used non-hydrostatic extension of shallow water theory is equivalent to a class of Boussinesq-type equations [7]. This involves augmenting the purely two-dimensional equation set by an additional equation for vertical velocity and an expression for the vertical (non-hydrostatic) pressure component, which yields additional terms in the momentum equation \((2)\) accounting for dispersion and vertical pressure profile effects.

In conclusion, while they are often considered outdated and oversimplified, shallow water equations still provide an interesting and rich field of research. I anticipate further new results and maybe even surprises from presentations at the upcoming SIAM Conference on Mathematical and Computational Issues in the Geosciences, to be held in Erlangen, Germany in September 2017.

¹ Funded by the Volkswagen Foundation

References
[1] Androsov, A., Behrens, J., & Danilov, S. (2011). Tsunami modelling with unstructured grids. Interaction between tides and tsunami waves. In E. Krause, Y. Shokin, M. Resch, D. Kröner, & N. Shokina (Eds.), Computational Science and High Performance Computing IV: Vol. 115. Notes on Numerical Fluid Mechanics and Multidisciplinary Design (pp. 191-206). Berlin: Springer.
[2] Behrens, J., & Dias, F. (2015). New computational methods in tsunami science, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 373.
[3] Behrens, J., & LeVeque, R.J. (2011, May). Modeling and simulating tsunamis with an eye to hazard mitigation. SIAM News, 44(4), p. 1.
[4] Dietrich, J.C., Tanaka, S., Westerink, J.J., Dawson, C.N., Luettich, R.A., Zijlema, M.,…Westerink, H.J. (2012). Performance of the unstructured-mesh, swan+adcirc model in computing hurricane waves and surge. Journal of Scientific Computing, 52, 468-497.
[5] Gruber, U., & Bartelt, P. (2007). Snow avalanche hazard modelling of large areas using shallow water numerical methods and GIS. Environmental Modelling & Software: Modelling, computer-assisted simulations, and mapping of dangerous phenomena for hazard assessment, 22(10), 1472-1481.
[6] Harig, S., Chaeroni, Pranowo, W.S., & Behrens, J. (2008). Tsunami simulations on several scales: Comparison of approaches with unstructured meshes and nested grids. Ocean Dynamics, 58, 429-440.
[7] Jeschke, A., Pedersen, G.K., Vater, S., & Behrens, J. (2017). Depth-averaged non-hydrostatic extension for shallow water equations with quadratic vertical pressure profile: Equivalence to boussinesq-type equations. Int. J. Numer. Meth. Fluids.
[8] LeVeque, R.J., & George, D.L. (2008). High-resolution finite volume methods for the shallow water equations with bathymetry and dry states. Advanced numerical models for simulating tsunami waves and runup, 10, 43-73.
[9] Özgen, I., Zhao, J., Liang, D., & Hinkelmann, R. (2016). Urban flood modeling using shallow water equations with depth-dependent anisotropic porosity. Journal of Hydrology: Part B, 541, 1165-1184.
[10] Popinet, S. (2012). Adaptive modelling of long-distance wave propagation and fine-scale flooding during the tohoku tsunami. Natural Hazards and Earth System Sciences, 12, 1213-1227.
[11] Saint-Venant, A. (1871). Théorie du mouvement non permanent des eaux, avec application aux crues des rivières et a l’introduction de marées dans leurs lits. Comptes Rendus des S ances de l Acad mie des Sciences, 73, 147, 154, & 237-240.
[12] Vater, S., Beisiegel, N., & Behrens, J. (2015). A limiter-based well-balanced discontinuous galerkin method for shallow-water flows with wetting and drying: One-dimensional case. Advances in Water Resources, 85, 1-13.
[13] Williamson, D.L., Drake, J.B., Hack, J.J., Jakob, R., & Swarztrauber, P.N. (1992). A standard test set for numerical approximations to the shallow water equations in spherical geometry. Journal of Computational Physics, 102, 211-224.

Jörn Behrens is a professor for numerical methods in geosciences at the University of Hamburg. He was trained as a mathematician and delivered his dissertation in adaptive atmospheric modeling at the University of Bremen, in Germany. After an assistant professorship at Technical University of Munich, and with an extended visit to the National Center for Atmospheric Research (NCAR) in Boulder, Colo., he took over the lead of the tsunami modeling group at Alfred-Wegener Institute in Bremerhaven, Germany. There he was responsible for developing the simulation component of the German-Indonesian Tsunami Early Warning System (GITEWS), now operational in Jakarta, Indonesia. His research interests cover adaptive numerical methods for multi-scale geoscientific applications, high-performance computing, and ocean and atmospheric modeling.