An energy and potential enstrophy conserving numerical scheme for the multi-layer shallow water equations with complete Coriolis force

https://doi.org/10.1016/j.jcp.2015.12.042Get rights and content

Highlights

  • Hamiltonian-derived scheme exactly conserves energy and potential enstrophy.

  • Energy/enstrophy conserved at machine accuracy for sufficiently small time step.

  • Can be used to incorporate full Coriolis force into existing isopycnal models.

Abstract

We present an energy- and potential enstrophy-conserving scheme for the non-traditional shallow water equations that include the complete Coriolis force and topography. These integral conservation properties follow from material conservation of potential vorticity in the continuous shallow water equations. The latter property cannot be preserved by a discretisation on a fixed Eulerian grid, but exact conservation of a discrete energy and a discrete potential enstrophy seems to be an effective substitute that prevents any distortion of the forward and inverse cascades in quasi-two dimensional turbulence through spurious sources and sinks of energy and potential enstrophy, and also increases the robustness of the scheme against nonlinear instabilities. We exploit the existing Arakawa–Lamb scheme for the traditional shallow water equations, reformulated by Salmon as a discretisation of the Hamiltonian and Poisson bracket for this system. The non-rotating, traditional, and our non-traditional shallow water equations all share the same continuous Hamiltonian structure and Poisson bracket, provided one distinguishes between the particle velocity and the canonical momentum per unit mass. We have determined a suitable discretisation of the non-traditional canonical momentum, which includes additional coupling between the layer thickness and velocity fields, and modified the discrete kinetic energy to suppress an internal symmetric computational instability that otherwise arises for multiple layers. The resulting scheme exhibits the expected second-order convergence under spatial grid refinement. We also show that the drifts in the discrete total energy and potential enstrophy due to temporal truncation error may be reduced to machine precision under suitable refinement of the timestep using the third-order Adams–Bashforth or fourth-order Runge–Kutta integration schemes.

Introduction

Geophysical fluid dynamics is concerned with the large-scale motions of stratified fluid on a rotating, near-spherical planet. In principle such motions are governed by the compressible Navier–Stokes equations as modified to include rotation, but simplified and approximate forms of these equations are widely used in both theoretical and computational models [1], [2]. At large scales viscous effects are relatively unimportant, and usually neglected. The resulting inviscid equations possess various conservation laws, for instance of energy, momentum, and potential vorticity [2]. The last property is perhaps the most important. Energy and momentum may be transported over large distances by waves, but potential vorticity is tied to fluid elements [3]. It is desirable for the simplified and approximated equations to share these conservation properties, and (as far as possible) for any numerical scheme to preserve them too.

Planetary scale fluid flow is typically dominated by the Coriolis force. However, almost all models of the ocean and atmosphere employ the so-called “traditional approximation” [4]. This retains only the contribution to the Coriolis force from the component of the planetary rotation vector that is locally normal to geopotential surfaces. The traditional approximation is typically justified on the basis that the vertical lengthscales of rotationally-dominated geophysical flows are much smaller than their horizontal lengthscales [5]. However, the “non-traditional” component of the Coriolis force is dynamically important in a wide range of geophysical flows [6]. In the ocean, nontraditional effects can substantially alter the structure and propagation properties of gyroscopic and internal waves [7], [8], [9], [10], [11], [12], the development and growth rate of inertial instabilities [13], [14], the instability of Ekman layers [15], high-latitude deep convection [16], the static stability of the water column [17], and the flow of abyssal waters close to the equator [18], [19], [20]. A parallel line of development has led to “deep atmosphere” models that include non-traditional effects [21], [22], [23], [24], [25], [26].

The focus of this article is a set of shallow water equations that include the full Coriolis force [27], [28]. Shallow water models approximate the ocean's stratification using a stack of layers of constant-density fluid. This is equivalent to a Lagrangian vertical discretisation of the water column using density as the vertical coordinate. Density or isopycnal coordinates are a natural choice for ocean models because mixing and transport across isopycnal surfaces in the ocean interior is weak, so large-scale flows are largely constrained to flow along density surfaces [29]. Models based on vertical coordinates aligned with surfaces of constant geopotential tend to produce spuriously large diabatic mixing due to the slight misalignment between grid cells and isopycnal surfaces [30]. Shallow water and isopycnal coordinate models are thus widely used to study the ocean circulation, both in conceptual modelling studies and comprehensive ocean models [31], [32], [33]. This isopycnal modelling framework has recently been extended to solve the non-hydrostatic three-dimensional stratified equations with the non-traditional components of the Coriolis force [34].

In this paper we derive an energy- and potential enstrophy-conserving numerical scheme for our extended shallow water equations that include the complete Coriolis force [27], [28]. There are three principal motivations for exactly preserving these conservation laws in our scheme:

  • 1.

    Ideally any numerical scheme should materially conserve potential vorticity, as this is a crucial dynamical variable in large-scale geophysical flows [3]. Material conservation of potential vorticity, Dq/Dt=0, implies conservation of all the integral moments of q,ddthc(q)dx=0, for any function c(q). These spatial integrals are much more amenable to approximation on a fixed grid. The most important is conservation of potential enstrophy, for which c(q)=q2. However, exact material conservation of potential vorticity is only possible via specialised Lagrangian schemes such as the particle-mesh method [35] or the contour-advective semi-Lagrangian method [36]. Thompson [37] showed that a vorticity evolution equation for an incompressible fluid that conserves energy and potential enstrophy, and satisfies reasonable properties such as spatial locality (for the streamfunction) and invariance under rotations, is inevitably of material conservation form. Although derived for partial differential equations in continuous space, this result suggests that conservation of energy and potential enstrophy on a fixed Eulerian grid provides a suitable substitute for material conservation of potential vorticity.

  • 2.

    Conservation of energy and potential enstrophy is crucial for the development of the simultaneous inverse cascade of energy and forward cascade of potential enstrophy that characterises quasi-two-dimensional turbulence [2]. Long-term integrations using non-conservative finite-difference schemes are prone to a spurious cascade of energy towards the scale of the grid spacing, where the spatial truncation error is largest [38]. Conserving the discrete analogue of total energy alone does not prevent the cascade of energy to small scales, and leads to an increase in potential enstrophy [39]. Thus long-time integrations require a scheme that conserves both energy and potential enstrophy, though even this is not sufficient to completely eliminate spurious spectral transfers of these quantities [40].

  • 3.

    Enforcing these conservation properties eliminates the possibility of nonlinear instabilities driven by spurious sources of energy or potential enstrophy at the grid scale. For example, Fornberg [41] showed that the only stable member of a class of semi-discrete centred finite difference schemes for the inviscid Burgers equation is the scheme that conserves a discrete energy.

Arakawa [38] constructed a finite difference scheme for the two-dimensional incompressible Euler equations, now known as the “Arakawa Jacobian”, that exactly conserves a discrete kinetic energy and a discrete enstrophy when formulated in continuous time. Arakawa and Lamb [42] later extended this approach to construct a finite difference scheme for the inviscid shallow water equations that exactly conserves total mass, energy, and potential enstrophy. The spurious energy cascades observed in non-conservative long-time integrations of the SWE (see point 2 above) are due to errors introduced in the spatial, rather than temporal, discretisation [39], as we shall confirm in §5. The Arakawa–Lamb scheme therefore discretises only the spatial derivatives of the traditional shallow water equations, yielding evolution equations in continuous time for the discrete particle velocities and layer thicknesses at spatial grid points. This semi-discrete finite difference scheme has since been extended to multiple layers [43], to fourth-order accuracy for the vorticity in the non-divergent limit [44], [45], [46], [47], to spherical geodesic grids [48], and to unstructured finite element meshes [49].

Arakawa and Lamb [42] began with a general finite difference discretisation of the shallow water equations on a staggered grid, known as the Arakawa C-grid, as sketched in Fig. 1. They wrote down a discretisation that contained many undetermined coefficients, including some unphysical Coriolis terms that were parallel, not perpendicular, to the fluid velocity. They obtained a scheme that conserves consistent spatially discrete approximations to the energy and potential enstrophy via direct algebraic manipulation of the discrete equations to determine suitable coefficients. In particular, the coefficients of the unphysical Coriolis terms become proportional to grid spacing squared. That a discrete scheme should exist that conserves not quadratic invariants like the Arakawa Jacobian [38], but a cubic quantity (energy) and a quotient (potential enstrophy) was long mysterious.

Salmon [47] interpreted the Arakawa–Lamb scheme as a particular discretisation of a non-canonical Hamiltonian formulation of the traditional shallow water equations, in which the evolution is expressed using a Hamiltonian functional and Poisson bracket. Symmetries of a Hamiltonian system correspond directly to conservation laws via Noether's theorem, making this a natural starting point for the development of conservative schemes. The equations of ideal fluid mechanics, including the shallow water equations, can be expressed as non-canonical Hamiltonian systems using space-fixed Eulerian variables [50], [2], [51], [52]. To keep the canonical form of Hamilton's equations one must either use Lagrangian variables that follow fluid elements, or introduce artificial Clebsch potentials.

The non-rotating, traditional, and our non-traditional shallow water equations all share the same Hamiltonian structure and Poisson bracket, provided one distinguishes between the particle velocity and the canonical momentum per unit mass [27], [28]. This distinction is often glossed over in the traditional shallow water equations because the two differ only by a function of space alone, so their time derivatives are identical. We use this insight to construct an analog of the Arakawa–Lamb scheme for our non-traditional shallow water equations by adapting Salmon's [47] Hamiltonian formulation. The main difficulty we encounter is that the layer thickness field enters the canonical momentum through the non-traditional terms, so we need to find a suitable average of the layer thickness field that is compatible with the placement of variables on the Arakawa C grid. Like the Arakawa Jacobian and Arakawa–Lamb schemes, we obtain a system of ordinary differential equations (ODEs) for the evolution in time of field values at grid points. Discrete analogs of the energy E and potential enstrophy P are exactly conserved under the exact evolution of this ODE system. To get a concrete numerical scheme we further discretise in time using the 4th order Runge–Kutta and 3rd order Adams–Bashforth schemes [53]. We investigate the time step that is necessary to reduce the temporal truncation errors in our discrete approximations of E and P down to machine precision.

Recently a similar Hamiltonian and Poisson bracket formulation has been employed to derive energy conserving schemes for the deep and shallow atmosphere equations [54], again using the canonical momenta to formulate the discrete Poisson bracket and ensure antisymmetry. However, these authors were unable to construct a discrete Poisson bracket that conserves potential enstrophy. A more general approach to constructing schemes that conserve energy and potential enstrophy formulates the shallow water equations as evolution equations for the potential vorticity, divergence, and layer thickness fields using discrete Nambu brackets [55], [56]. This formulation requires an elliptic inversion to compute the velocity field, which is computationally expensive but allows for a straightforward incorporation of boundary conditions.

The structure of this paper is as follows. In §2 we briefly review the shallow water equations with complete Coriolis force and their Hamiltonian formulation [27], [28]. In §3 we review Salmon's Hamiltonian derivation of the Arakawa–Lamb scheme for the traditional shallow water equations (§3.1), and extend it to include the complete Coriolis force in a single layer (§3.2) and in multiple layers (§3.3). In §4 we present a series of numerical verification experiments, and in §5 we demonstrate that energy and potential enstrophy can be conserved to machine precision if the time-step is sufficiently short. Finally, in §6 we summarise our results and provide concluding remarks.

Section snippets

Shallow water equations with complete Coriolis force

The derivation of our numerical scheme begins with the non-canonical Hamiltonian formulation of the traditional shallow water equations, henceforth abbreviated as SWE, as described in [50], [51], [52]. The single and multi-layer non-traditional shallow water equations (henceforth NTSWE) were recently derived in [27], [28]. For clarity we first present the non-canonical Hamiltonian formulation of the single-layer NTSWE.

We begin by writing the single layer NTSWE [27], [28] in the vector invariant

Spatial discretisation

In this section we derive an energy- and potential enstrophy-conserving spatial discretisation for the NTSWE. We begin with a brief review of the Arakawa–Lamb scheme [42] for the single-layer SWE and its Hamiltonian formulation [47], then extend the scheme to the NTSWE, and to multiple layers.

Numerical experiments

We first nondimensionalise our equations by writingx=Rdxˆ,u(k)=cuˆ(k),h(k)=Hhˆ(k),hb=Hhˆb,(Ω(x),Ω(y),Ω(z))=Ω(Ωˆ(x),Ωˆ(y),Ωˆ(z)), where a hat aˆ denotes a dimensionless variable. We use the layer thickness scale H, gravitational acceleration g, and planetary rotation rate Ω to construct the gravity wave speed c, Rossby deformation radius Rd, and non-traditional parameter δ [18], [22], [27]c=gH,Rd=c2Ω,δ=H/Rd=2ΩH/c. The last of these is the aspect ratio based on the deformation radius and

Temporal truncation error

The convergence studies described in §4 demonstrate that the numerical scheme derived in §3 achieves the expected second-order accuracy under spatial grid refinement. However, the exact conservation of the discrete total energy E and potential enstrophy P only holds in the semi-discrete formulation in which the values of u, v, h at grid points evolve according to an ODE system in continuous time. In this section we show that the changes E and P due to the finite accuracy of the time integration

Conclusion

We have extended the energy and potential enstrophy conserving properties of the Arakawa–Lamb [42] scheme to the multi-layer non-traditional shallow water equations (NTSWE) that include a depth-averaged treatment of the complete Coriolis force [27], [28]. These conservation properties help to prevent a spurious forward energy cascade to the finest resolved scales, and suppress nonlinear numerical instabilities. They seem essential for correctly simulating the inverse energy cascade in

Acknowledgements

Some of the computations employed the Advanced Research Computing (ARC) facility at the University of Oxford [70]. The authors thank the three anonymous reviewers for their constructive comments on the manuscript.

References (70)

  • R. Salmon

    Lectures on Geophysical Fluid Dynamics

    (1998)
  • P. Muller

    Ertel's potential vorticity theorem in physical oceanography

    Rev. Geophys.

    (1995)
  • C. Eckart

    Variation principles of hydrodynamics

    Phys. Fluids

    (1960)
  • N.A. Phillips

    The equations of motion for a shallow rotating atmosphere and the traditional approximation

    J. Atmos. Sci.

    (1966)
  • T. Gerkema et al.

    Geophysical and astrophysical fluid dynamics beyond the traditional approximation

    Rev. Geophys.

    (2008)
  • H. van Haren et al.

    Gyroscopic waves in the Mediterranean Sea

    Geophys. Res. Lett.

    (2005)
  • H. van Haren

    Sharp near-equatorial transitions in inertial motions and deep-ocean step-formation

    Geophys. Res. Lett.

    (2005)
  • T. Gerkema et al.

    Near-inertial waves in the ocean: beyond the “traditional approximation”

    J. Fluid Mech.

    (2005)
  • T. Gerkema et al.

    Near-inertial waves on the non-traditional β-plane

    J. Geophys. Res.

    (2005)
  • T. Gerkema et al.

    Non-traditional reflection of internal waves from a sloping bottom, and the likelihood of critical reflection

    Geophys. Res. Lett.

    (2006)
  • A.L. Stewart et al.

    Multilayer shallow water equations with complete Coriolis force. Part 2. Linear plane waves

    J. Fluid Mech.

    (2012)
  • B.L. Hua et al.

    Inertial nonlinear equilibration of equatorial flows

    J. Fluid Mech.

    (1997)
  • A. Colin de Verdière

    The stability of short symmetric internal waves on sloping fronts: beyond the traditional approximation

    J. Phys. Oceanogr.

    (2012)
  • S. Leibovich et al.

    The influence of the horizontal component of Earth's angular velocity on the instability of the Ekman layer

    J. Fluid Mech.

    (1985)
  • J. Marshall et al.

    Open-ocean convection: observations, theory, and models

    Rev. Geophys.

    (1999)
  • A.L. Stewart et al.

    Multilayer shallow water equations with complete Coriolis force. Part 3. Hyperbolicity and stability under shear

    J. Fluid Mech.

    (2013)
  • A.L. Stewart et al.

    Cross-equatorial channel flow with zero potential vorticity under the complete Coriolis force

    IMA J. Appl. Math.

    (2012)
  • R. Müller

    A note on the relation between the “traditional approximation” and the metric of the primitive equations

    Tellus A

    (1989)
  • A.A. White et al.

    Dynamically consistent, quasi-hydrostatic equations for global models with a complete representation of the Coriolis force

    Q. J. R. Meteorol. Soc.

    (1995)
  • A.A. White

    A view of the equations of meteorological dynamics and various approximations

  • A.A. White et al.

    Consistent approximate models of the global atmosphere: shallow, deep, hydrostatic, quasi-hydrostatic and non-hydrostatic

    Q. J. R. Meteorol. Soc.

    (2005)
  • M. Tort et al.

    Dynamically consistent shallow-atmosphere equations with a complete Coriolis force

    Q. J. R. Meteorol. Soc.

    (2014)
  • M. Tort et al.

    Usual approximations to the equations of atmospheric motion: a variational perspective

    J. Atmos. Sci.

    (2014)
  • P.J. Dellar et al.

    Shallow water equations with a complete Coriolis force and topography

    Phys. Fluids

    (2005)
  • A.L. Stewart et al.

    Multilayer shallow water equations with complete Coriolis force. Part 1. Derivation on a non-traditional beta-plane

    J. Fluid Mech.

    (2010)
  • Cited by (30)

    • A mass conservative, well balanced, tangency preserving and energy decaying method for the shallow water equations on a sphere

      2022, Journal of Computational Physics
      Citation Excerpt :

      Li et al. [19] proposed a finite difference weighted essentially non-oscillatory (WENO) interpolation-based scheme to maintain well balancedness (the exact C-property) as well as to mitigate the Gibbs oscillations. Energy conserving semidiscrete spatial discretizations for the inviscid shallow water equations were developed based on mixed finite-element methods [20], finite element exterior calculus on the sphere [10], Hamiltonian formulation [27,29,34], and mixed mimetic spectral element [16,17]. Hoang et al. [13] developed a fully discrete scheme which can preserve mass conservation and potential vorticity by combining the C-grid staggering method [25] with the strong stability preserving Runge–Kutta (SSP-RK) method.

    • Parallel exponential time differencing methods for geophysical flow simulations

      2021, Computer Methods in Applied Mechanics and Engineering
    View all citing articles on Scopus
    View full text