# A simple and transparent method for improving the energetics and thermodynamics of seawater approximations: Static energy asymptotics (SEA)

Rémi Tailleux, Thomas Dubos

► To cite this version:

Rémi Tailleux, Thomas Dubos. A simple and transparent method for improving the energetics and thermodynamics of seawater approximations: Static energy asymptotics (SEA). Ocean Modelling, 2024, 188, 5, pp. 439-461. 10.1016/j.ocemod.2024.102339. insu-04847484

HAL Id: insu-04847484

https://insu.hal.science/insu-04847484v1

Submitted on 19 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A simple and transparent method for improving the energetics and thermodynamics of seawater approximations: Static energy asymptotics (SEA)

Rémi Tailleux<sup>a,*</sup>, Thomas Dubos<sup>b</sup>

## ARTICLE INFO

**Keywords:**

Sound-proof approximations

Seawater Boussinesq/Anelastic approximation

Static energy

Energetically and thermodynamically

consistent models

Thermodynamic potentials

## ABSTRACT

The static energy encodes all possible information about the thermodynamics and potential energy (and all related forces) of stratified geophysical fluids. In this paper, we develop a systematic methodology, called static energy asymptotics, that exploits this property for constructing energetically and thermodynamically consistent sound-proof approximations of the equations of motion. By approximating the static energy to various orders of accuracy, two main families of approximations are (re-)derived and discussed: the pseudo-incompressible (PI) approximation and the anelastic (AN) approximation. For all approximations, the background and available potential energies (in Lorenz sense) can be constructed to match their exact counterparts as closely as feasible and to be expressible in terms of the exact (as opposed to ad-hoc) thermodynamic potentials. For hydrostatic motions, the AN approximation (of which the Boussinesq approximation is a special case) has the same structure as that of legacy Seawater Boussinesq primitive equations. The energetics of such models could therefore be made transparently traceable to that of the full Navier-Stokes equations at little to no additional cost, thus allowing them to take full advantage of the Gibbs Sea Water (GSW) library developed as part of the new thermodynamic standard for seawater TEOS-10.

## 1. Introduction

Since their inception in the late 60's, numerical ocean models have almost exclusively relied on the so-called Seawater Boussinesq approximation (SBA thereafter). Unlike the standard Boussinesq approximation, the SBA makes use of the full nonlinear equation of state for seawater and therefore retains the adiabatic and isohaline compressibility effects resulting from its pressure dependence. Doing so is essential because these nonlinearities appear to be crucial for correctly simulating the relative layering of ocean water masses, e.g., Nycander et al. (2015), IOC, SCOR, IAPSO (2010). The nonlinearities most important dynamically are: the *cabbeling nonlinearity*, which is associated with the temperature dependence of the thermal expansion coefficient and responsible for the densification upon mixing, and the *thermobaric nonlinearity*, which is associated with the pressure dependence of the thermal expansion coefficient and responsible for making colder parcels more compressible than warmer ones (McDougall, 1987).

Although the SBA has received much attention, its accuracy, exact justification, as well as its energetic and thermodynamic consistency, have remained poorly understood and sources of confusion. Conceptually, the construction and justification of the SBA used in numerical

ocean models a priori involves a two-step process: (1) the first one pertaining to the coarse-graining of the un-averaged equations of motion; (2) the second one pertaining to the sound-proofing of the equations. So far, most approaches have generally assumed that (2) should be carried out before (1), as it is easiest to implement and understand, but the opposite view also exists. For instance, McDougall et al. (2002) have argued that the divergence-free velocity field used by the SBA should be interpreted as the Favre averaged velocity field $\vec{v}^\rho = \overline{\rho\nabla\cdot\vec{v}}$, thus regarding the SBA as being primarily the result of the coarse-graining procedure rather than of the Boussinesq sound-proofing step. Clarifying the issue is important for establishing the accuracy of the SBA, because in McDougall's interpretation, the continuity equation $\nabla \cdot \vec{v}^\rho = 0$ is exact in a steady-state, whereas it is only approximate in the conventional interpretation.

In this paper, our main aim is to clarify the key ingredients (and those that are not) of any construction of energetically and thermodynamically consistent approximations, defined here as approximations that admit an energy conservation principle and thermodynamic potentials traceable in a rigorous and physically transparent way to

<figure><img src="image_3.png" /></figure>

<footnote><sup>a</sup> Department of Meteorology, University of Reading, Whiteknights road, Earley Gate, Reading, RG6 6ET, Berkshire, United Kingdom</footnote>

<footnote><sup>b</sup> IPSL, Lab. Meteorologie Dynamique, Ecole Polytechnique, Palaiseau, France</footnote>

<footnote>* Corresponding author.
E-mail address: r.g.j.tailleux@reading.ac.uk (R. Tailleux).</footnote>
their exact counterparts, and whose description of irreversible processes/diabatic effects remains in agreement with the second law of thermodynamics. To that end, we address the issue in the context of the un-averaged compressible Navier-Stokes equations, as how coarse-graining affects energetics and thermodynamics is currently much less understood. Note that in the literature, the concept of energetically and thermodynamically consistent models is also occasionally understood as models predicting the evolution of coarse-grained motions whose turbulent closures are consistent with energy conservation and the second law of thermodynamics. Those models may also include prognostic equations for subgrid scale energy reservoirs, e.g., Eden et al. (2014), Eden (2015, 2016).

Following relatively recent progress (Young, 2010; Tailleux, 2010b, 2012; Eden et al., 2014; Eden, 2015, 2016), the SBA is now understood to admit a well defined energy conservation principle, although one that requires the construction of ad-hoc thermodynamic potentials. One of the main results of this paper is that this undesirable feature stems from the un-necessary character of some of the approximations/simplifications made by the SBA. Indeed, we find that it is always possible to construct SBA-like approximations whose thermodynamics can be expressed in terms of the exact thermodynamic potentials. The approximations that we propose, unlike the legacy SBA, can therefore take full advantage of the Gibbs Seawater Library (GSW) developed as part of the new international thermodynamic standard TEOS-10 (IOC, SCOR, IAPSO, 2010).

It has long been observed that the SBA can be regarded as being isomorphic to the fully compressible hydrostatic equations of motion written in pressure coordinates (de Szoeke and Samelson, 2002). This result is important, because it implies that SBA-model can in principle be made fully compressible and therefore more accurate simply by re-interpreting height as pressure and implementing a few other changes without fundamentally altering the architecture of existing codes, as discussed and tested by Losch et al. (2004) for instance. The use of pressure coordinates is much less natural for the oceans than it is for the atmosphere, however, since it requires considering a prognostic equation for the bottom pressure (which is difficult to observe) rather than for the free surface (observable by satellite altimetry). For this reason, there is a strong incentive to continue searching for alternative ways to improve on the SBA that retain standard coordinates. To that end, there appears to be a wealth of sound-proof approximations to choose from. The anelastic system (Ogura and Phillips, 1962; Lipps and Hemler, 1982) has been initially derived for dry air modelled as an ideal perfect gas and adiabatic conditions. Scinocca and Shepherd (1992) have obtained a Hamiltonian formulation of the anelastic system, which has later been extended to accommodate arbitrary multicomponent thermodynamics (Pauluis, 2008). Vasil et al. (2013), Tort and Dubos (2014b), Cotter and Holm (2014) have developed variational formulations for these systems that enable to trace back their conservation properties to symmetries of an adequate Lagrangian. Durran (1989, 2008) has proposed a different approximation, known as the pseudo-incompressible (PI) approximation in an atmospheric context. Arbitrary thermodynamics and diabatic effects (i.e. molecular conduction and diffusion) have been consistently introduced in this system (Klein and Pauluis, 2012). In an oceanic context, Dewar et al. (2016) have proposed two types of “semi-compressible” approximations with a consistent treatment of diabatic effects, the type-I approximation coinciding with the PI system and type-II with the anelastic system.

Most recently, Eldred and Gay-Balmaz (2021) have devised variational formulations of the PI and anelastic systems that extend Vasil et al. (2013) and Tort and Dubos (2014b) by including diabatic effects that previously had to be dealt with directly at the level of the equations of motion (Pauluis, 2008; Klein and Pauluis, 2012; Dewar et al., 2016). Based on the existing literature, the problem of how to construct energetically and thermodynamically consistent sound-proof approximations to the equations of motion appears to be largely understood, at least in the context of the un-averaged equations of motion. In

most cases, however, such approaches require a significant amount of mathematical background, e.g. differential geometry and variational principles with non-holonomic constraints. In the context of the SBA, the methods developed by Tailleux (2012), Eden et al. (2014), Eden (2015, 2016) to achieve energetics and thermodynamic consistency are arguably much simpler, but unfortunately not sufficiently general to be easily extended to the pseudo-incompressible or anelastic approximations. As a result, there is still a strong incentive to look for alternative approaches that can make the process of constructing consistent approximations more physically intuitive while keeping the maths as elementary as feasible. The main aim of this paper is to reformulate the mathematical apparatus of Eldred and Gay-Balmaz (2021) in a much simpler unifying formalism, where the main vehicle of the various levels of approximations is the so-called static energy. This concept, which should be more familiar to both oceanographers and atmospheric scientists, is found to allow for a transparent discussion of both the energetics and accuracy of Boussinesq-like approximations. The hope is that by demystifying energetics and thermodynamic consistency issues it might facilitate the adoption of more straightforwardly consistent and/or more accurate approximations by numerical ocean modellers and others.

The paper is organised as follows: Section 2 introduces the topic by first reminding the reader of where to locate all the information about the thermodynamics and energetics in the compressible Navier-Stokes equations. This leads us to identify the static energy as the most natural variable in which to encode all thermodynamic information about the fluid, and therefore the one to approximate. Section 3 shows how to construct the pseudo-incompressible approximation. Section 4 discusses the construction of a modernised anelastic approximation. Section 5 discusses the construction of modernised versions of the Boussinesq and anelastic approximations for seawater within a single unifying framework. Section 6 discusses the impact of the approximation on the partitioning of the potential energy into available potential energy (APE) and background potential energy (BPE), as per Tailleux (2018) local formulation of Lorenz (1955) theory of available potential energy. Section 7 summarises and discusses the results.

## 2. Energy conservation and structure of static energy

To understand how to construct energetically and thermodynamically consistent approximations, it is first important to understand how the equations of motion encode the information about thermodynamics and potential energy. It is also important to understand how this information might get partially lost or scrambled in the simplification or approximation process underlying the construction of idealised models or of sound-proof approximations, or even through the use of non-canonical variables. Indeed, this is necessary for identifying possible information recovery strategies. The following sections discuss and illustrate these issues.

### 2.1. Remark on energy conservation in the Navier-Stokes equations

We take as our ground-truth and reference model the Navier-Stokes equations for compressible seawater approximated as a two-constituent fluid, which may be written

$$ \frac{D\nu}{Dt} + 2\Omega \times \nu + \frac{1}{\rho}\nabla p = -\nabla\Phi + \mathbf{F}, \qquad (1) $$

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0 \quad \text{or equivalently} \quad \nabla \cdot \mathbf{v} = \frac{D}{Dt} \ln \nu, \qquad (2) $$

$$ \frac{D\eta}{Dt} = \dot{\eta} = -\frac{1}{\rho}\nabla \cdot (\rho \mathbf{J}_{\eta}) + \dot{\eta}_{irr} \qquad (3) $$

$$ \frac{DS}{Dt} = \dot{S} = -\frac{1}{\rho}\nabla \cdot (\rho \mathbf{J}_S), \qquad (4) $$

$$ T = T(\eta, S, p) = \frac{\partial h}{\partial \eta}, \quad \mu = \mu(\eta, S, p) = \frac{\partial h}{\partial S} \qquad (5) $$
$$v = v(\eta, S, p) = \frac{\partial h}{\partial p} \qquad (6)$$

Here, $h = h(\eta, S, p)$ is the specific enthalpy, $v = (u, v, w)$ is the three-dimensional velocity field, $\rho$ is density, $v = 1/\rho$ is the specific volume, $p$ is pressure, $\Phi = gz$ is the geopotential, $\Omega$ is Earth rotation vector, $g$ is the acceleration of gravity, $z$ is height, $S$ is salinity, $\eta$ is specific entropy, with $J_S$ and $J_\eta$ the molecular fluxes of salt and entropy, respectively. The term $\dot{\eta}_{irr} \ge 0$ represents the irreversible rate of specific entropy production, which the second law of thermodynamics requires to be non-negative. Finally, $\mathbf{F} = \rho^{-1}\nabla \cdot \mathbf{S}$ is the viscous force, where $\mathbf{S}$ is the deviatoric stress tensor.

Consistency with the second law of thermodynamics is associated with the evolution equations for specific entropy and salinity (3) and (4) through the specification of the phenomenological laws for the molecular diffusive fluxes $J_\eta$ and $J_S$, as well as with the constraint that $\dot{\eta}_{irr} \ge 0$ be positive. From non-equilibrium thermodynamics theory, we know that the molecular diffusive fluxes act to relax the fluid towards a resting state ($v = 0$) with homogeneous in-situ temperature $T = \text{constant}$ and relative chemical potential $\mu = \text{constant}$. The simplest way to achieve this is via assuming that the molecular fluxes of $\eta$ and $S$ be linearly related to the generalised forces $\nabla T/T$ and $\nabla \mu/T$ as follows

$$J_{\eta} = -L_{\eta\eta} \frac{\nabla T}{T} - L_{\eta s} \frac{\nabla \mu}{T} \qquad (7)$$

$$J_S = -L_{\eta s} \frac{\nabla T}{T} - L_{ss} \frac{\nabla \mu}{T} \qquad (8)$$

with $L_{\eta\eta}$, $L_{ss}$, $L_{\eta s}$ and $L_{s\eta}$ to be determined empirically from laboratory measurements. Energy conservation can be shown to impose that the following quantity

$$\rho[\mathbf{v} \cdot \mathbf{F} + T\dot{\eta} + \mu\dot{S}] = -\nabla \cdot (\rho J_E) \qquad (9)$$

be equal to the divergence of an appropriately defined energy flux $J_E$. This in turn constrains the non-viscous irreversible entropy production to take the form

$$\dot{\eta}_{irr}^{\text{diff}} = L_{\eta\eta} \frac{|\nabla T|^2}{T^2} + L_{ss} \frac{|\nabla \mu|^2}{T^2} + (L_{s\eta} + L_{\eta s}) \frac{\nabla T \cdot \nabla \mu}{T^2} \qquad (10)$$

which is non-negative provided that

$$L_{\eta\eta} > 0, \quad L_{ss} > 0, \quad L_{\eta\eta} L_{ss} - L_{\eta s} L_{s\eta} \ge 0. \qquad (11)$$

Finally, the Onsager reciprocity relationship allows for the simplification $L_{\eta s} = L_{s\eta}$. As shown in the following, the main effect of a sound-proof approximation is to modify the form of the in-situ temperature and relative chemical potential entering the phenomenological laws (7) and (8) but not the laws themselves.

As is well known, in addition to being consistent with the second law of thermodynamics, the Navier-Stokes equations are built to satisfy the law of energy conservation, the conserved total energy per unit mass being $E_{tot} = E_k + E_p$, with $E_k = v^2/2$ the kinetic energy (KE), and

$$E_p = \Phi(z) + h - \frac{p}{\rho} = \Phi(z) + h - p \frac{\partial h}{\partial p} \qquad (12)$$

the potential energy (PE), itself the sum of gravitational potential energy $\Phi(z) = gz$ and internal energy $e = h - p/\rho$. Note that the Navier-Stokes equations contain only information about the three partial derivatives of $h$,

$$\frac{\partial h}{\partial p} = v = \frac{1}{\rho}, \quad \frac{\partial h}{\partial \eta} = T, \quad \frac{\partial h}{\partial S} = \mu, \qquad (13)$$

meaning that they are only sufficient to recover $h$ up to an irrelevant integration constant. (Note though that while $v$ and $T$ are unambiguously defined, $\mu$ is only defined up to an arbitrary constant; it follows that $h$ is fundamentally defined only up to a linear function of $S$, see IOC, SCOR, IAPSO, 2010). Partial loss of information about $h$ commonly occurs in many systems of equations used in practice, however. This is the case, for instance, when restricting attention to adiabatic and isohaline

motions ($\dot{\eta} = \dot{S} = 0$), for which knowledge of the thermohaline derivatives of $h$ ($T$ and $\mu$) is no longer formally needed for integrating the equations forward in time. In this case, the Navier-Stokes equations reduce to the Euler equations (assuming that viscous processes are also neglected), and only contain enough information to recover $h$ (or $e$) up to an indeterminate function of $(\eta, S)$. In the literature, the latter are often referred to as energy Casimirs, e.g., Shepherd (1993).

Physically, it is important to remark that specific enthalpy acts as a thermodynamic potential encoding all possible thermodynamic information only if formulated in terms of its canonical variables $(\eta, S, p)$ (Alberty, 1994). It follows that thermodynamic information may also be lost when using expressions of specific enthalpy formulated in terms of non-canonical variables such as potential temperature $\theta$ or Conservative Temperature $\Theta$ in place of entropy as $h = \dot{h}(S, \theta, p)$ or $\dot{h}(S, \Theta, p)$, as is usually preferred in oceanographic practice. To avoid losing information in that case, it becomes necessary to also supply the passage relations $\eta = \dot{\eta}(S, \theta, p)$ or $\eta = \dot{\eta}(S, \Theta, p)$ allowing one to reconstruct $h$ in terms of its canonical variables, e.g., IOC, SCOR, IAPSO (2010). This is needed, for instance, to compute variables such as in-situ temperature $T$, whose expression is

$$T = \frac{\partial h}{\partial \eta} = \frac{\partial \dot{h}}{\partial \theta} / \frac{\partial \dot{\eta}}{\partial \theta}, \qquad (14)$$

which clearly requires knowledge of both $\dot{h}$ and $\dot{\eta}$. In practice, the loss of thermodynamic information affecting widely used systems of equations such as the standard Seawater Boussinesq approximation used by ocean modellers and discussed next may therefore have multiple origins that one needs to be aware of, as it is obviously a pre-requisite for identifying information recovery strategies.

## 2.2. Motivating problem: seawater Boussinesq approximation

To make the above ideas more concrete, it is useful to discuss the particular case of the Seawater Boussinesq approximation (SBA) used in a majority of numerical Ocean General Circulation Models (OGCMs) and the main focus of this paper. Its governing equations are

$$\frac{\mathrm{D}\mathbf{u}}{\mathrm{D}t} + f \mathbf{k} \times \mathbf{u} + \nabla_h \left( \frac{\delta p}{\rho_b} \right) = \mathbf{F}_h \qquad (15)$$

$$\frac{\partial}{\partial z} \left( \frac{\delta p}{\rho_b} \right) = b_{\text{bou}}, \qquad (16)$$

$$b_{\text{bou}} = b_{\text{bou}}(S, \theta, z) = -\frac{g(\rho - \rho_b)}{\rho_b} \qquad (17)$$

$$\nabla \cdot (\rho_b \mathbf{v}) = 0 \qquad (18)$$

$$\frac{\mathrm{D}S}{\mathrm{D}t} = -\nabla \cdot \mathbf{F}_S, \quad \frac{\mathrm{D}p}{\mathrm{D}t} = -\nabla \cdot \mathbf{F}_\theta \qquad (19)$$

$$\rho = \rho(S, \theta, p_b(z)) \qquad (20)$$

where $\rho_b$ is the constant Boussinesq reference density, which in the NEMO OGCM is generally chosen to be $\rho_b = 1026 \text{ kgm}^{-3}$, $p_b = -\rho_b gz$ is the Boussinesq pressure, $\delta p = p - p_b(z)$. As $\rho_b$ is constant, the continuity equation is of course equivalent to the usual form $\nabla \cdot \mathbf{v} = 0$. Remaining notations are: $\mathbf{u} = (u, v)$ is the horizontal component of the coarse-grained 3D velocity field $\mathbf{v} = (u, v, w)$; $\nabla_h$ is the horizontal component of the $\nabla$ operator; $f$ is the Coriolis parameter; $\mathbf{k}$ is the unit normal vector pointing upward; $\mathbf{F}_h$ is the horizontal component of the turbulent viscous force; $\mathbf{F}_S$ and $\mathbf{F}_\theta$ are the parameterised turbulent fluxes of salt and heat.

Although (15)–(20) has formed the basis for a majority of OGCMs, that they admit an energy conservation principle was only clarified relatively recently by Young (2010), Tailleux (2012), and Eden (2015) among others, building upon previous ideas by Ingersoll (2005) and Pauluis (2008). Starting from the kinetic energy equation,

$$\rho_b \frac{\mathrm{D} \mathbf{u}^2}{\mathrm{D} t} + \nabla \cdot (\delta p \mathbf{v}) = \rho_b b_{\text{bou}}(S, \theta, z) \frac{\mathrm{D} z}{\mathrm{D} t} + \rho_b \mathbf{F}_h \cdot \mathbf{u}, \qquad (21)$$
obtained by taking the scalar product of (15) with **u** combined with (16) multiplied by Dz/Dt, the key step here is to recognise that the term $b_{bou}Dz/Dt$ is the Lagrangian derivative of a pseudo potential energy $E_{p,bou}$, so that (21) may be rewritten in the form

$$ \rho_b \frac{D}{Dt} \left( \frac{\mathbf{u}^2}{2} + E_{p,bou} \right) + \nabla \cdot (\delta p \mathbf{v}) = \rho_b \left( Y_S \frac{DS}{Dt} + Y_\theta \frac{D\theta}{Dt} + \mathbf{F}_h \cdot \mathbf{u} \right) \quad (22) $$

where the Boussinesq pseudo potential energy $E_{p,bou}$ and the thermodynamic efficiencies $Y_S$ and $Y_\theta$ are defined by

$$ E_{p,bou} = - \int_{z_*}^{z} b_{bou}(S, \theta, \bar{z}) \, d\bar{z}, \quad (23) $$

$$ Y_S = - \int_{z_*}^{z} \frac{\partial b_{bou}}{\partial S}(S, \theta, \bar{z}) \, d\bar{z}, \quad (24) $$

$$ Y_{\theta} = - \int_{z_*}^{z} \frac{\partial b_{bou}}{\partial \theta}(S, \theta, \bar{z}) \, d\bar{z}, \quad (25) $$

where $z_*$ is an arbitrary reference depth. Eq. (22) defines an energy conservation principle for $\mathbf{u}^2/2 + E_{p,bou}$ for inviscid motions in the absence of diabatic sources/sinks of $\theta$ and $S$, but is unsatisfactory because $E_{p,bou}$ represents only one chunk of the total potential energy, whose determination is rendered ambiguous by the arbitrariness of $z_*$ and whose link to its exact counterpart is unclear. While progress has been made towards resolving these issues, existing ideas still remain largely ad-hoc and lacking in generality, however, thus motivating the considerations developed in subsequent sections.

## 2.3. Importance of static energy

To ensure that the full energy conservation can be recovered from the equations of motion, we find it essential to explicitly articulate the connections between aspects of the momentum balance equations and potential energy. In this paper, this is achieved by expressing the pressure-geopotential gradient force $\rho^{-1}\nabla p + \nabla\Phi$, potential energy $E_p$, and dynamics/thermodynamics coupling in terms of the static energy

$$ \Sigma_{FC}(\eta, S, p, \Phi) = h(\eta, S, p) + \Phi, \quad (26) $$

by means of the following relations:

$$ \frac{1}{\rho} \nabla p + \nabla \Phi = \frac{\partial \Sigma_{FC}}{\partial p} \nabla p + \frac{\partial \Sigma_{FC}}{\partial \Phi} \nabla \Phi = \mathbf{F}_{dyn}, \quad (27) $$

$$ E_p = \Phi + h - \frac{p}{\rho} = \Sigma_{FC} - p \frac{\partial \Sigma_{FC}}{\partial p}, \quad (28) $$

$$ \underbrace{\nabla \Sigma_{FC} - \frac{\partial \Sigma_{FC}}{\partial \eta} \nabla \eta - \frac{\partial \Sigma_{FC}}{\partial S} \nabla S}_{\text{thermodynamics}} = \underbrace{\frac{\partial \Sigma_{FC}}{\partial p} \nabla p + \frac{\partial \Sigma_{FC}}{\partial \Phi} \nabla \Phi}_{\text{dynamics}}, \quad (29) $$

where in the case of a fully compressible fluid, the partial derivatives of $\Sigma_{FC}$ are

$$ \frac{\partial \Sigma_{FC}}{\partial p} = v = \frac{1}{\rho}, \quad \frac{\partial \Sigma_{FC}}{\partial \Phi} = 1, \quad \frac{\partial \Sigma_{FC}}{\partial \eta} = T, \quad \frac{\partial \Sigma_{FC}}{\partial S} = \mu. \quad (30) $$

Notice that in Eq. (26), $(p, \eta, S, \Phi)$ are independent coordinates in an augmented state space that extends the standard thermodynamic state space $(p, \eta, S)$. By a standard abuse of notations, we use the same notation for these independent variables and for the corresponding fields.

In the literature, the identity (29) is sometimes called the Crocco-Vazsonyi theorem (Crocco, 1937; Vazsonyi, 1945) (referred as the CVT thereafter) and is key for connecting the dynamics (27) to the thermodynamics. Physically, it implies that the momentum balance equations, whose standard dynamical form is

$$ \frac{Dv}{Dt} + 2\Omega \times v + \frac{\partial \Sigma_{FC}}{\partial p} \nabla p + \frac{\partial \Sigma_{FC}}{\partial \Phi} \nabla \Phi = \mathbf{F} \quad (31) $$

may also be written in thermodynamic form as follows:

$$ \frac{Dv}{Dt} + 2\Omega \times v + \nabla \Sigma_{FC} = \mathbf{F} + \frac{\partial \Sigma_{FC}}{\partial \eta} \nabla \eta + \frac{\partial \Sigma_{FC}}{\partial S} \nabla S, \quad (32) $$

referred to as the Crocco equations thereafter. As shown throughout this paper, the three identities (27), (28), and (29) greatly facilitate understanding of energy conservation issues. The key idea of this paper is to argue that to understand the energetics and thermodynamic consistency of any given approximation to the equations of motion, it is essential to understand how said approximation affects these identities.

To proceed, let us first establish that energy conservation is actually attached to the structure of $\Sigma_{FC}$ rather than to its explicit form, as this will prove essential to obtain a general principle applicable to both the approximated and non-approximated equations of motions. In the following, we therefore only assume that: (1) $\Sigma_{FC}$ is some general function of $(\eta, S, p, \Phi)$; (2) that the density entering the continuity equation Eq. (2) relates to the pressure derivative of $\Sigma_{FC}$ via

$$ \rho \frac{\partial \Sigma_{FC}}{\partial p} = 1. \quad (33) $$

Thus, taking the scalar product of (32) with $\mathbf{v}$, using the facts that by definition $\mathbf{v} \cdot \nabla C = DC/Dt - \partial C/\partial t = \dot{C} - \partial C/\partial t$ for any scalar field $C(x, t)$ and that $\partial\Phi/\partial t = 0$, yields

$$ \begin{aligned} \frac{D}{Dt} \left( \frac{\mathbf{v}^2}{2} + \Sigma_{FC} \right) &= \mathbf{F} \cdot \mathbf{v} + \frac{\partial \Sigma_{FC}}{\partial \eta} \dot{\eta} + \frac{\partial \Sigma_{FC}}{\partial S} \dot{S} + \frac{\partial \Sigma_{FC}}{\partial t} \dot{t} - \frac{\partial \Sigma_{FC}}{\partial \eta} \frac{\partial \eta}{\partial t} - \frac{\partial \Sigma_{FC}}{\partial S} \frac{\partial S}{\partial t} \\ &= \mathbf{F} \cdot \mathbf{v} + \frac{\partial \Sigma_{FC}}{\partial \eta} \dot{\eta} + \frac{\partial \Sigma_{FC}}{\partial S} \dot{S} + \frac{\partial \Sigma_{FC}}{\partial p} \frac{\partial p}{\partial t}. \end{aligned} \quad (34) $$

Next, multiplying the result by $(\partial\Sigma_{FC}/\partial p)^{-1} = \rho$, accounting for the continuity equation Eq. (2), yields

$$ \frac{\partial}{\partial t} (\rho E_{total}) + \nabla \cdot \left[ \rho \left( \frac{\mathbf{v}^2}{2} + \Sigma_{FC} \right) \mathbf{v} \right] = \rho \left( \mathbf{F} \cdot \mathbf{v} + \frac{\partial \Sigma_{FC}}{\partial \eta} \dot{\eta} + \frac{\partial \Sigma_{FC}}{\partial S} \dot{S} \right) \quad (35) $$

with $E_{total} = E_k + E_p$, with $E_k = \mathbf{v}^2/2$ and

$$ E_p = \Sigma_{FC} - p \frac{\partial \Sigma_{FC}}{\partial p}. \quad (36) $$

Eq. (35) makes it clear that $E_{total}$ is conserved in the absence of viscous ($\mathbf{F} = 0$) and diabatic ($\dot{\eta} = \dot{S} = 0$) effects. For total energy conservation to hold in the general case, it is necessary to make the extra assumption that the right-hand side of (35) may be written as the divergence of some flux $J_E$ including viscous and diffusive effects, viz.,

$$ \rho \left( \mathbf{F} \cdot \mathbf{v} + \frac{\partial \Sigma_{FC}}{\partial \eta} \dot{\eta} + \frac{\partial \Sigma_{FC}}{\partial S} \dot{S} \right) = - \nabla \cdot (\rho J_E). \quad (37) $$

Physically, Eq. (37) in turn constrains the form of $\dot{\eta}_{lrr}$ in terms of the molecular fluxes $J_\eta$ and $J_s$ as well as with the gradients of $\nabla T$ and $\nabla \mu$ to ensure consistency with the second law of thermodynamics, e.g., see Pauluis (2008), Tailleux (2010a) and Woods (1975). Note that (35) may alternatively be written in the form

$$ \frac{\partial(\rho B_{FC})}{\partial t} + \nabla \cdot (\rho B_{FC} \mathbf{v}) = \rho \frac{DB_{FC}}{Dt} = -\nabla \cdot (\rho J_E) + \frac{\partial p}{\partial t} \quad (38) $$

where

$$ B_{FC} = \frac{\mathbf{v}^2}{2} + \Sigma_{FC} \quad (39) $$

is the standard Bernoulli function attached to $\Sigma_{FC}$. Eq. (38) is the appropriate form of the energy conservation principle for establishing the Bernoulli theorem, namely that $B_{FC}$ is conserved along adiabatic, isohaline, inviscid, steady fluid parcel trajectories.

## 3. Pseudo-Incompressible (PI) approximation

To summarise, we established that under the conditions stated above, the equations of motion formulated in terms of the static energy $\Sigma_{FC}(\eta, S, p, \Phi)$ as per (31) or (32), together with the continuity equation Eq. (2), the entropy and salinity budgets (3)–(4) and the constraint (33) conserve the total energy $E_{total} = E_k + E_p$, with the potential energy being given by

$$ E_p = \Sigma_{FC} - p \frac{\partial \Sigma_{FC}}{\partial p} \quad (40) $$
regardless of the exact form of $\Sigma_{FC}$. As a consequence, the same result must hold for any approximation to the equations of motion obtained by any structure-preserving approximation of $\Sigma_{FC}$, thus providing a systematic procedure to derive such approximations, baptised here static energy asymptotics. In this section, this approach is illustrated by re-deriving the pseudo-incompressible (PI) approximation (Durran, 1989, 2008; Dewar et al., 2016; Eldred et al., 2022) via replacing $\Sigma_{FC}$ by an approximation $\Sigma_{PI}$ defined in terms of the leading order terms of a Taylor series expansion around the reference pressure $p_R(\Phi)$ as follows:

$$ \Sigma_{FC} = \underbrace{h(\eta, S, p_R(\Phi)) + v(\eta, S, p_R(\Phi))(p - p_R(\Phi)) + \Phi + R_{PI}}_{\Sigma_{PI}}, \quad (41) $$

where $\delta p = p - p_R(\Phi)$ is the assumed small pressure anomaly and $R_{PI} = O(\delta p^2)$ the residual term. Importantly, Eq. (41) shows that $\Sigma_{PI} = \Sigma_{PI}(\eta, S, p, \Phi)$ retains the same functional form as $\Sigma_{FC}$ and hence that the above results must apply. Evaluating the $(\eta, S, p, \Phi)$ derivatives of $\Sigma_{PI}$ yields

$$ \frac{\partial \Sigma_{PI}}{\partial p} = v(\eta, S, p_R(\Phi)) = v_* = \frac{1}{\rho_*}, \quad (42) $$

$$ \frac{\partial \Sigma_{PI}}{\partial \Phi} = 1 + \frac{\rho_R(\Phi)\delta p}{\rho_*^2 c_{s*}^2}, \quad (43) $$

$$ \frac{\partial \Sigma_{PI}}{\partial \eta} = T(\eta, S, p_R(\Phi)) + \delta p \frac{\partial v_*}{\partial \eta} = T_* + \delta p \frac{\partial v_*}{\partial \eta} = T_{PI}, \quad (44) $$

$$ \frac{\partial \Sigma_{PI}}{\partial S} = \mu(\eta, S, p_R(\Phi)) + \delta p \frac{\partial v_*}{\partial S} = \mu_* + \delta p \frac{\partial v_*}{\partial S} = \mu_{PI}, \quad (45) $$

and show that these remain close to the derivatives of $\Sigma_{FC}$, with relatively minor impact on $\rho$, $T$, and $\mu$, the main effect being the actual pressure $p$ being replaced by the reference pressure $p_R(\Phi)$ and the presence of a small correction $O(\delta p)$, where $\rho_R(\Phi) = -d p_R/d\Phi$ denotes the reference density in hydrostatic equilibrium with $p_R$, while $c_{s*} = (\partial \rho/\partial p)^{-1/2}(\eta, S, p_R(\Phi))$ is the speed of sound evaluated at $p_R(\Phi)$. Generally, the star subscript refers to quantities evaluated at the reference pressure $p_R(\Phi)$ instead of the full pressure $p$. Note that the linearity of $\Sigma_{PI}$ in $p$ makes $v_* = \partial \Sigma_{PI}/\partial p$ independent of $p$, which filters sound waves out due to the resulting decoupling of density and pressure. As a result, the above approximation affects the Crocco theorem, momentum balance equation, and continuity equation as follows:

$$ \nabla \Sigma_{PI} - T_{PI} \nabla \eta - \mu_{PI} \nabla S = \frac{1}{\rho_*} \nabla p + \frac{\partial \Sigma_{PI}}{\partial \Phi} \nabla \Phi, \quad (46) $$

$$ \frac{Dv}{Dt} + 2\Omega \times v + \frac{1}{\rho_*} \nabla p + \frac{\partial \Sigma_{PI}}{\partial \Phi} \nabla \Phi = \mathbf{F}, \quad (47) $$

$$ \frac{\partial \rho_*}{\partial t} + \nabla \cdot (\rho_* \mathbf{v}) = 0, \quad (48) $$

noting that, per (33) applied to $\Sigma_{PI}$ instead of $\Sigma_{FC}$, $\rho$ must also be replaced by $\rho_*$ in the tracer equations for specific entropy and salinity (3) and (4). Since $\Sigma_{PI}$ has the same functional structure as $\Sigma_{FC}$, it follows that the results derived in the previous section apply and hence that the associated approximations of motion conserve the total energy $E_{total} = E_k + E_{p,PI}$, with

$$ E_{p,PI} = \Sigma_{PI} - \frac{p}{\rho_*} = h(\eta, S, p_R) - \frac{p_R}{\rho_*} + \Phi = e(\eta, S, p_R) + \Phi \quad (49) $$

provided that the diabatic terms $\dot{\eta}$ and $\dot{S}$ are constrained to satisfy

$$ \rho_* [\mathbf{F} \cdot \mathbf{v} + T_{PI} \dot{\eta} + \mu_{PI} \dot{S}] = -\nabla \cdot (\rho_* \mathbf{J}_E). \quad (50) $$

In that case, (49) shows that the approximation to the potential energy $E_{p,PI}$ is close to its exact counterpart, and hence that the energetics of the PI approximation is traceable to that of the fully compressible equations. Eq. (49) shows that for the PI approximation to be consistent with the second law, it is sufficient to replace $T$, $\mu$, and $\rho$ by $T_{PI}$, $\mu_{PI}$, and $\rho_*$ in the expressions for the diabatic terms $\dot{\eta}$ and $\dot{S}$, as previously found by Klein and Pauluis (2012) and Eldred and Gay-Balmaz (2021).

### 3.1. Remarks on the choice of $p_R(\Phi)$

Physically, the accuracy of the PI approximation is determined by the magnitude of the residual $R_{PI}$ in (41), which may be rewritten in the following equivalent forms:

$$ \begin{aligned} R_{PI} &= h(\eta, S, p) - h(\eta, S, p_R(\Phi)) - v(\eta, S, p_R(\Phi))(p - p_R(\Phi)) \\ &= \int_{p_R(\Phi)}^{p} [v(\eta, S, \bar{p}) - v(\eta, S, p_R(\Phi))] d\bar{p} = \int_{p_R(\Phi)}^{p} \int_{p_R(\Phi)}^{\bar{p}} v_p(\eta, S, \bar{p}) d\bar{p} d\bar{p}. \end{aligned} \quad (51) $$

Because $v_p = -1/(\rho_*^2 c_{s*}^2) < 0$ where $c_s$ is the speed of sound, it follows that $R_{PI}$ is also negative. Moreover, it is also easily shown that $R_{PI}$ is quadratic in $\delta p$ at leading order, viz.,

$$ R_{PI} \approx -\frac{(p - p_R(\Phi))^2}{2\rho_*^2 c_{s*}^2}, \quad (52) $$

and hence that its magnitude is directly controlled by the choice of $p_R(\Phi)$ and its distance to $p$. The reader familiar with the local theory of APE will may recognise that $R_{PI}$ is approximately equal to the opposite of the available compressible energy (ACE) denoted by $\Pi_1$ in Tailleux (2018). As a consequence,

$$ \Sigma_{PI} = h(\eta, S, p_R(\Phi)) + \Phi + \frac{\delta p}{\rho_*} = \Sigma_{FC} - R_{PI} \geq \Sigma_{FC}, \quad (53) $$

which establishes, perhaps counter-intuitively, that $\Sigma_{PI}$ represents an overestimate of the true $\Sigma_{FC}$. The construction of $\Sigma_{PI}$ requires identifying a suitable reference pressure $p_R(\Phi)$, but how best to do so in practice is rarely discussed. The fact that $R_{PI}$ is sign definite suggests, however, that $p_R(\Phi)$ should be defined to minimise the volume integral of $R_{PI}$. In the case where the denominator $\rho_*^2 c_{s*}^2$ can be treated as approximately constant, this would amount to defining $p_R(\Phi)$ as the horizontal-mean pressure.

### 3.2. Degenerate energy conservation principles and uniqueness issues

As is well known, the energy associated with the law of energy conservation for fully compressible fluids described by the NSE is $E_{tot} = E_k + E_p$, and therefore represents the fundamental form of energy against which to assess the energy conservation principles attached to any approximation to the equations of motion. The main aim of this section is to demonstrate that a degenerate energy conservation principle may occasionally be obtained in the context of the approximations discussed in this paper, which is important to be understood and be aware of for correctly interpreting some results of the literature.

As shown below, an alternative energy conservation principle may be shown to exist in the PI approximation because

$$ \Sigma_{PI} = h(\eta, S, p_R(\Phi)) + v(\eta, S, p_R(\Phi))\delta p + \Phi = \hat{\Sigma}_{PI}(\eta, S, \delta p, \Phi) \quad (54) $$

(from Eq. (41)) can alternatively be interpreted as a function of $(\eta, S, \delta p, \Phi)$ (denoted with a hat). Since $\delta p = p - p_R(\Phi)$ is a function of $p$ and $\Phi$ only, the pressure-geopotential gradient force can therefore be written equivalently as

$$ \frac{\partial \Sigma_{PI}}{\partial p} \nabla p + \frac{\partial \Sigma_{PI}}{\partial \Phi} \nabla \Phi = \frac{\partial \hat{\Sigma}_{PI}}{\partial \delta p} \nabla \delta p + \frac{\partial \hat{\Sigma}_{PI}}{\partial \Phi} \nabla \Phi, \quad (55) $$

the $(\delta p, \Phi)$ partial derivatives of $\hat{\Sigma}_{PI}$ being given by

$$ \frac{\partial \hat{\Sigma}_{PI}}{\partial \delta p} = \frac{\partial \Sigma_{PI}}{\partial p} = \frac{1}{\rho_*}, \quad (56) $$

$$ \frac{\partial \hat{\Sigma}_{PI}}{\partial \Phi} = 1 - \frac{\rho_R(\Phi)}{\rho_*} + \frac{\delta p \rho_R(\Phi)}{\rho_*^2 c_{s*}^2}, \quad (57) $$

while in contrast its thermodynamic derivatives remain unaffected

$$ \frac{\partial \hat{\Sigma}_{PI}}{\partial \eta} = \frac{\partial \Sigma_{PI}}{\partial \eta} = T_{PI}, \quad \frac{\partial \hat{\Sigma}_{PI}}{\partial S} = \frac{\partial \Sigma_{PI}}{\partial S} = \mu_{PI}. \quad (58) $$
From (55), it follows that the PI momentum balance equations may alternatively be written in the form

$$ \frac{D\mathbf{v}}{Dt} + 2\Omega \times \mathbf{v} + \frac{1}{\rho_{\star}} \nabla \delta p + \frac{\partial \hat{\Sigma}_{PI}}{\partial \Phi} \nabla \Phi = \mathbf{F}, \qquad (59) $$

which from the same structural arguments developed previously implies the existence of a degenerate form of energy conservation for the pseudo total energy $E_{\text{total}} = E_k + E_{\text{p,pseudo}}$, with

$$ E_{\text{p,pseudo}} = \hat{\Sigma}_{PI} - \delta p \frac{\partial \hat{\Sigma}}{\partial p} = \hat{\Sigma}_{PI} - \frac{\delta p}{\rho_{\star}} = h(\eta, S, p_R(\Phi)) + \Phi. \qquad (60) $$

However,

$$ \Sigma_{PI} - p \frac{\partial \Sigma_{PI}}{\partial p} \neq \hat{\Sigma}_{PI} - \delta p \frac{\partial \hat{\Sigma}_{PI}}{\partial p} \qquad (61) $$

so that $E_{\text{p,pseudo}}$ and $E_{\text{p,PI}}$ differ, specifically by the term $p_R(\Phi)/\rho_{\star}$. However, as the latter quantity may be shown to satisfy the conservation law

$$ \rho_{\star} \frac{D}{Dt} \left( \frac{p_R(\Phi)}{\rho_{\star}} \right) = \nabla \cdot (p_R(\Phi)\mathbf{v}), \qquad (62) $$

it follows that the conservation principle for $E_k+E_{\text{p,pseudo}}$ simply results from adding (62) to the conservation principle for the correct energy $E_k+E_{\text{p,PI}}$. This particular example highlights that achieving a flux form conservation law for an energy-like quantity is not sufficient to establish the energy conservation principle satisfied by the approximation considered: one also needs to be able to relate the energy-like quantity obtained to that of the fully-compressible energy, a point that many authors do not appear to be aware of, see Eq. (26) of Dewar et al. (2016) for instance.

## 4. Anelastic approximation (AN) revisited

In the PI approximation, the density $\rho_{\star} = \rho(\eta, S, p_R(\Phi))$<fnref n="1" /> becomes decoupled from pressure but still responds to changes in salinity and entropy. It therefore remains time-dependent, unlike in the Boussinesq approximation. In the anelastic (AN) approximation, which we revisit in this section, the density is further approximated to be a static function of $\Phi$ only, $\rho_{\star}(\Phi)$. The attendant static energy $\Sigma_{AN}$ is obtained by further approximating $\Sigma_{PI}$ as follows

$$ \Sigma_{AN} = h(\eta, S, p_R(\Phi)) + \frac{\delta p}{\rho_{\star}} + \Phi, \qquad (63) $$

and is obtained by replacing $\rho_{\star} = \rho(\eta, S, p_R(\Phi))$ by $\rho_{\star}(\Phi)$ in (41). Per (33), the continuity equation applies to $\rho_{\star}(\Phi)$, viz.: $\nabla \cdot (\rho_{\star}v) = 0$. Physically, this approach unifies (i) the standard anelastic approximation where $\rho_{\star} = \rho_R(\Phi)$ is a general function of $\Phi$ (Ogura and Phillips, 1962; Pauluis, 2008; Eldred and Gay-Balmaz, 2021) and (ii) $\rho_{\star} = \rho_R = \text{constant}$, corresponding to a modernised form of Boussinesq approximation.

The AN approximation $\Sigma_{AN}$ retains the same $(\eta, S, p, \Phi)$ functional dependence characterising $\Sigma_{FC}$ and $\Sigma_{PI}$. Its partial derivatives are:

$$ \frac{\partial \Sigma_{AN}}{\partial p} = \frac{1}{\rho_{\star}}, \qquad (64) $$

$$ \frac{\partial \Sigma_{AN}}{\partial \Phi} = 1 + \rho_R (v_{\star} - v_{\star}) + \delta p \frac{dv_{\star}}{d\Phi}, \qquad (65) $$

$$ \frac{\partial \Sigma_{AN}}{\partial \eta} = T_{\star}, \qquad (66) $$

$$ \frac{\partial \Sigma_{AN}}{\partial S} = \mu_{\star}, \qquad (67) $$

(with $v_{\star} = \rho_{\star}^{-1}$) and can be seen to remain close to the derivatives of $\Sigma_{PI}$, albeit with a loss of accuracy due to the absence of the small correction $O(\delta p)$ and to $\rho_{\star} \neq \rho_{\star}^-$.

The fact that $\Sigma_{AN}$ has the same structure as $\Sigma_{FC}$ or $\Sigma_{PI}$ means that provided that its diabatic and viscous terms satisfy the constraint

$$ \rho_{\star} \left[ \mathbf{F} \cdot \mathbf{v} + T_{\star} \dot{\eta} + \mu_{\star} \dot{S} \right] = - \nabla \cdot (\rho_{\star} \mathbf{J}_{E}) \qquad (68) $$

analogous to (9), it must conserve the total energy $E_k+E_{\text{p,AN}}$, with the AN potential energy being given by

$$ E_{\text{p,AN}} = \Sigma_{AN} - p \frac{\partial \Sigma_{AN}}{\partial p} = h(\eta, S, p_R(\Phi)) - \frac{p_R}{\rho_{\star}} + \Phi = E_{\text{p,PI}} + p_R(v_{\star} - v_{\star}). \qquad (69) $$

Eq. (69) shows that $E_{\text{p,AN}}$ differs only from $E_{\text{p,PI}}$ by the energy error term $\Delta E_{\text{p,spurious}} = p_R(v_{\star} - v_{\star})$. This error can be minimised by defining $v_{\star}$ as the temporally and horizontally averaged value of $v_{\star}$ for instance.

Eqs. (64)–(67) show that the CVT can be written in the standard form

$$ \nabla \Sigma_{AN} - T_{\star} \nabla \eta - \mu_{\star} \nabla \mu = \frac{1}{\rho_{\star}} \nabla p + \frac{\partial \Sigma_{AN}}{\partial \Phi} \nabla \Phi. \qquad (70) $$

However, a more convenient and concise expression of the CVT can be obtained by taking advantage of the possibility of regarding $\Sigma_{AN} = \hat{\Sigma}_{AN}(\eta, S, \delta p/\rho_{\star}, \Phi)$ as a function of $(\eta, S, \delta p/\rho_{\star}, \Phi)$ instead, for which the $(\eta, S, \Phi)$ derivatives can be shown to be

$$ \frac{\partial \hat{\Sigma}_{AN}}{\partial \eta} = \frac{\partial \Sigma_{AN}}{\partial \eta} = T(\eta, S, p_R(\Phi)) = T_{\star}, \qquad (71) $$

$$ \frac{\partial \hat{\Sigma}_{AN}}{\partial S} = \frac{\partial \Sigma_{AN}}{\partial S} = \mu(\eta, S, p_R(\Phi)) = \mu_{\star}, \qquad (72) $$

$$ \frac{\partial \hat{\Sigma}_{AN}}{\partial \Phi} = 1 - \rho_R(\Phi)v(\eta, S, p_R(\Phi)) = 1 - \frac{\rho_R(\Phi)}{\rho_{\star}} = -\frac{b_{AN}}{g} \qquad (73) $$

where $b_{AN}$ is recognised as the ‘buoyancy’ defined relative to $\rho_R(\Phi)$, which differs from the traditional Boussinesq buoyancy as further discussed below. It follows that the CVT (70) may be rewritten in the form

$$ \nabla \hat{\Sigma}_{AN} - T_{\star} \nabla \eta - \mu_{\star} \nabla S = \nabla \left( \frac{\delta p}{\rho_{\star}} \right) - \frac{b_{AN}}{g} \nabla \Phi. \qquad (74) $$

Eq. (74) differs from (70) in that the density $\rho_{\star}$ is now absorbed within the anomalous pressure gradient term, while only the buoyancy $b_{AN}$ multiplies $\nabla\Phi$, thus making it possible to rewrite the standard form of momentum balance in the form

$$ \frac{D\mathbf{v}}{Dt} + 2\Omega \times \mathbf{v} + \nabla \left( \frac{\delta p}{\rho_{\star}} \right) - \frac{b_{AN}}{g} \nabla \Phi = \mathbf{F}. \qquad (75) $$

The resulting diabatic anelastic system coincides with the one derived by Eldred and Gay-Balmaz (2021), Eq. (4.36). A key advantage is to make the form of hydrostatic balance similar to that of the standard Boussinesq approximation even in the case where $\rho_{\star}$ depends on $z$, as further discussed in the next section.

## 5. Modernised anelastic and Boussinesq seawater approximations

As established previously, the AN approximation provides an energetically and thermodynamically consistent model with realistic thermodynamic potentials that closely resemble their exact counterparts. In the context of hydrostatic ocean modelling that has been traditionally studied using legacy SBA primitive equations models discussed in Section 2(b), the AN approximation takes the form (reverting to using $z$ rather than $\Phi$, which is more standard in oceanographic practice)

$$ \frac{D\mathbf{u}}{Dt} + f\mathbf{k} \times \mathbf{u} + \nabla_h \left( \frac{\delta p}{\rho_{\star}} \right) = \mathbf{F}_h, \qquad (76) $$

$$ \frac{\partial}{\partial z} \left( \frac{\delta p}{\rho_{\star}} \right) = b_{AN}, \qquad (77) $$

<footnote n="1">¹ By density, we always mean the density field that obeys the continuity equation Eq. (2).</footnote>
$$b_{AN} = -g(1 - \rho_R v_\star), \qquad (78)$$

$$\nabla \cdot (\rho_\frac{3}{2} v) = 0, \qquad (79)$$

$$\frac{DS}{Dt} = F_S, \qquad \frac{DS}{Dt} = F_\theta, \qquad (80)$$

where $\rho_\frac{3}{2}(z)$ does not need to coincide with $\rho_R(z) = -p'_R(z)/g$. In practice, however, it seems natural to equate the two to avoid using multiple reference densities. While in the standard SBA the reference pressure has been commonly approximated as $p_R(z) = -\rho_b g z$, with $\rho_b$ = constant, arguments developed in next section suggest that $p_R(z)$ and $\rho_R(z)$ could be advantageously defined in terms of the reference profiles entering the local theory of available potential energy (APE).

Mathematically, the AN model (76)–(80) can be verified to have essentially the same structure as that of the legacy SBA model (15)–(20), and hence that it can be solved algorithmically and numerically in exactly the same way regardless of whether $\rho_\frac{3}{2}$ is taken as constant or function of $z$. The subtle difference between (76)–(80) and (15)–(20) is that in (78) $b_{AN}$ represents the non-approximated buoyancy advocated by Young (2012) and Tailleux (2012). Further advantages of the AN model are: (1) the form of potential energy entering its energy conservation principle, viz.,

$$E_{p,AN} = h(S, \theta, p_R(z)) + \Phi(z) - \frac{p_R(z)}{\rho_R(z)} \qquad (81)$$

is traceable to its exact counterpart without the need to introduce ad-hoc thermodynamic potentials, an undesirable feature of the legacy SBA models; (2) the physical basis for specifying the reference density profile $\rho_R(z)$ that it uses is much clearer than that for specifying the constant Boussinesq reference density $\rho_b$ used by legacy SBA models. For this reason, the above AN model is more capable of exploiting the new capabilities offered by the Gibbs Sea Water library developed as part of TEOS-10 to its full extent. Although $\rho_\frac{3}{2}$ could also be chosen as a constant reference Boussinesq density $\rho_b$, this simplification offers no advantage over using a depth-dependent $\rho_R(z)$ consistent with $p_R(z)$. This clearly shows, therefore, that some of the approximations made as part of the legacy Seawater Boussinesq approximation deteriorate its energetics without leading to real simplifications or computational benefits compared to the modernised version proposed here and to the anelastic system.

## 6. Dynamics/thermodynamics partitioning of energy

It is now well established that energy comes in at least two different flavours, heat-like versus work-like, or available versus non-available, or available (APE) and background (BPE), as justified by Lorenz (1955) theory of available potential energy (APE) for instance. Since according to thermodynamics, ‘work’ and ‘heat’ represent fundamentally distinct forms of energy, which in APE theory are regarded as dynamically active and inert respectively, it follows that any discussion of energetics and thermodynamics consistency cannot be complete without understanding how sound-proof approximations individually affect the active and passive forms of energy.

To express the above idea in the present framework, we assume that it is possible to meaningfully decompose the static energy into dynamically active ($\Sigma_{dyn}$) and passive ($\Sigma_{heat}$) components as follows:

$$\Sigma_{FC}(\eta, S, p, \Phi) = \Sigma_{dyn}(\eta, S, p, \Phi) + \Sigma_{heat}(\eta, S). \qquad (82)$$

Since $\Sigma_{heat}$ is independent of $p$ and $\Phi$, it follows that only $\Sigma_{dyn}$ is dynamically relevant for predicting the pressure-geopotential gradient force entering the standard form of momentum balance equations,

$$\frac{1}{\rho} \nabla p + \nabla \Phi = \frac{\partial \Sigma_{dyn}}{\partial p} \nabla p + \frac{\partial \Sigma_{dyn}}{\partial \Phi} \nabla \Phi, \qquad (83)$$

thus justifying the dynamically irrelevant character of $\Sigma_{heat}$. In the following, we briefly discuss the two main approaches for defining

$\Sigma_{dyn}$ and $\Sigma_{heat}$ proposed so far, and how these are affected by the most drastic AN approximation. Eq. (83) shows that to be physically acceptable, any construction of $\Sigma_{dyn}$ should only require knowledge of $\partial\Sigma/\partial p$ and $\partial\Sigma/\partial\Phi$ (regardless of how $\Sigma$ is defined), possibly combined with knowledge of the whole stratification at some point in time (such as provided by initial conditions).

For the sake of clarity and consistency with existing literature, we switch back to using height/depth $z$ as vertical coordinate.

### 6.1. Available potential energy (APE) theory

From a theoretical viewpoint, the most rigorous and physically-based partitioning of $\Sigma_{FC}$ is arguably the one rooted in APE theory. In this approach, $\Sigma_{dyn}$ is defined as the departure of $\Sigma_{FC}$ from its equilibrium value defined as the one it would have in a suitably defined notional reference state characterised by the reference pressure $p_0(z)$ and $\rho_0(z)$, that is,

$$\Sigma_{heat} = h(\eta, S, p_0(z_R)) + g z_R \qquad (84)$$

where $z_R = z_R(\eta, S)$ is the fluid parcel's reference position. Eq. (84) thus represents a local form of background potential energy (BPE), which in Lorenz (1955) was only defined as a volume-integrated quantity.

Because in a state of rest, the density of a fluid parcel must match that of the background reference state, $z_R$ must be a solution of the so-called level of neutral buoyancy (LNB) equation (Tailleux, 2013, 2018)

$$\rho(\eta, S, p_0(z_R)) = \rho_0(z_R). \qquad (85)$$

With $\Sigma_{heat}$ defined as per (84), we have

$$\Sigma_{dyn} = \Sigma_{FC} - \Sigma_{heat} = \Pi_1 + \Pi_2 + \frac{p - p_0(z)}{\rho} \qquad (86)$$

where $\Pi_1$ and $\Pi_2$ are so-called available compressible energy (ACE) density and APE density respectively,

$$\Pi_1 = h(\eta, S, p) - h(\eta, S, p_0(z)) + \frac{p_0(z) - p}{\rho} \qquad (87)$$

$$\Pi_2 = h(\eta, S, p_0(z)) - h(\eta, S, p_0(z_R)) + g(z - z_R). \qquad (88)$$

As shown by Tailleux (2018), these may be written as

$$\Pi_1 = \int_{p_0(z)}^{p} \int_{p}^{p'} \frac{\partial v}{\partial p}(\eta, S, p'') dp'' dp' \qquad (89)$$

$$\Pi_2 = - \int_{z_R}^{z} b(\eta, S, z') dz' \qquad (90)$$

where $b(\eta, S, z) = -g[1 - \rho_0(z)\nu(\eta, S, p_0(z))]$ is the buoyancy defined relative to the reference density $\rho_0(z)$. This establishes, therefore, that $\Pi_1$, $\Pi_2$, and therefore $\Sigma_{dyn}$ can be constructed only from the knowledge of $\nu(\eta, S, p)$, $g$, as well as from one state of the system from which $p_0(z)$ and $\rho_0(z)$ are constructed (for instance by using the adiabatic and isohaline re-arrangement technique discussed by Saenz et al., 2015).

Let us now seek to understand how the AN approximation affects the dynamics/thermodynamics partitioning of $\Sigma_{AN}$,

$$\Sigma_{AN} = h(\eta, S, p_R(z)) + g z + \frac{\delta p}{\rho_\frac{3}{2}}. \qquad (91)$$

By considering the value of $\Sigma_{AN}$ in the Lorenz background reference state, in which $z = z_R$ and $p = p_0(z_R)$, it is easily realised that $\Sigma_{AN, heat}$ can be made to coincide with its exact counterpart (84) simply by choosing $p_R(z) = p_0(z)$, thus providing a natural solution to the problem of how to define $p_R(z)$ in the PI and AN approximations. With such a choice, the dynamical component of $\Sigma_{AN}$ becomes

$$\Sigma_{AN,dyn} = \Sigma_{AN} - \Sigma_{AN,heat} = \Pi_2 + \frac{p - p_0(z)}{\rho_\frac{3}{2}}, \qquad (92)$$

where the expression for $\Pi_2$ is the same as its exact counterpart. By comparing (92) with its exact counterpart (86), it is seen that the AN
approximation results in the disappearance of the ACE component $\Pi_1$, while also approximating the compressible work term $(p - p_0(z))/\rho$ by $(p - p_0(z))/\rho_\ddagger$. To avoid the proliferation of reference density profiles, the most logical choice is to use $\rho_\ddagger = \rho_0(z)$, which provides a natural solution for how to choose $\rho_\ddagger$.

## 6.2. Potential/dynamic enthalpy partitioning

We briefly remark that another kind of partitioning is based on defining $\Sigma_{\text{heat}}$ as the 'potential' value of $\Sigma_{\text{FC}}$ referenced to the ocean surface at $z = 0$, viz.,

$$ \Sigma_{\text{heat}} = h(\eta, S, p_R(0)) + \Phi(0), \qquad (93) $$

(de Szoek, 2000; Young, 2010; Nycander, 2010) where $p_R(z)$ is an arbitrary hydrostatic reference pressure field defining the reference density field via $\rho_R(z) = -p'_R(z)/g$. If $p_R(0) = p_a$ and $\Phi(0) = 0$, where $p_a$ is the mean surface atmospheric pressure, $\Sigma_{\text{heat}}$ then coincides with Mc-Dougall (2003)'s potential enthalpy, the currently recommended variable for defining heat content in the oceans following the adoption of TEOS-10 (IOC, SCOR, IAPSO, 2010; Pawlowicz et al., 2012). With $\Sigma_{\text{heat}}$ defined as per (93), the expression for $\Sigma_{\text{dyn}}$ can easily be verified to be

$$ \Sigma_{\text{dyn}} = \Pi_1 + h^{\ddagger}(\eta, S, z) + \frac{p - p_R(z)}{\rho} \qquad (94) $$

with $\Pi_1 = h(\eta, S, p) - h(\eta, S, p_R(z)) + (p_R(z) - p)/\rho \ge 0$ the ACE defined in terms of $p_R(z)$ rather than $p_0(z)$, and $h^{\ddagger}$ the so-called dynamic enthalpy, defined by

$$ h^{\ddagger} = h(\eta, S, p_R(z)) - h(\eta, S, p_R(0)) + \Phi(z) - \Phi(0) = - \int_0^z b(\eta, S, z') dz' \quad (95) $$

with $b(\eta, S, z) = -g[1 - \rho_R(z)v(\eta, S, p_R(z))]$ the buoyancy defined relative to the reference density $\rho_R(z)$. It is easily seen that (93) and (94) only differ from their APE-based counterparts by the use of a constant reference level located at the surface in place of the parcels' resting position, as well as in the choice of reference density/pressure profiles. Similarly as for the APE-based BPE, the AN approximation does not affect the thermodynamic component of $\Sigma_{AN}$, so that $\Sigma_{AN,\text{heat}} = \Sigma_{\text{heat}}$. The AN approximation only affects the dynamical component by getting rid of $\Pi_1$ in (94), so that

$$ \Sigma_{AN,dyn} = h^{\ddagger} + \frac{\delta p}{\rho_{\ddagger}}. \qquad (96) $$

The potential/dynamic enthalpy based partitioning of $\Sigma_{AN}$, however, lacks any of the advantageous attributes and clear physical interpretation of the BPE/APE based partitioning; its usefulness is therefore questionable.

## 7. Summary and discussion

The main result of this paper is that it is possible to write the equations of motion for a general compressible two-constituent stratified fluid, as well as for a large class of thermodynamically and energetically consistent sound-proof approximations of such equations, under the generic form

$$ \Sigma = \Sigma(\eta, S, p, \Phi), \qquad (97) $$

$$ B = \frac{\partial \Sigma}{\partial \Phi}, \quad \rho^{-1} = v = \frac{\partial \Sigma}{\partial p}, \quad T = \frac{\partial \Sigma}{\partial \eta}, \quad \mu = \frac{\partial \Sigma}{\partial S}, \qquad (98) $$

$$ \frac{Dv}{Dt} + 2\Omega \times v = -v \nabla p + B \nabla \Phi + \mathbf{F}, \qquad (99) $$

$$ \nabla \cdot \mathbf{v} = \frac{D}{Dt} \ln v, \qquad (100) $$

$$ \frac{D\eta}{Dt} = -v \nabla \cdot (\rho \mathbf{J}_{\eta}) + \dot{\eta}_{lrr} = \dot{\eta}, \qquad (101) $$

$$ \frac{DS}{Dt} = -v \nabla \cdot (\rho \mathbf{J}_s) = \dot{S}, \qquad (102) $$

$$ \mathbf{F} \cdot \mathbf{v} + T\dot{\eta} + \mu\dot{S} = -v\nabla \cdot (\rho\mathbf{J}_E), \qquad (103) $$

$$ T\mathbf{J}_{\eta} = -L_{\eta\eta}\nabla T - L_{\eta s}\nabla\mu, \qquad (104) $$

$$ T\mathbf{J}_s = -L_{\eta s}\nabla T - L_{ss}\nabla\mu. \qquad (105) $$

As shown in this paper, such equations conserve total mass defined as the volume integral of $(\partial\Sigma/\partial p)^{-1}$, total salinity defined as the volume integral of $(\partial\Sigma/\partial p)^{-1} S$ and total energy defined as the volume integral of $(\partial\Sigma/\partial p)^{-1} E_{tot}$, where $E_{tot} = v^2/2 + \Sigma - \rho\partial\Sigma/\partial p$. Furthermore, they are consistent with the second law of thermodynamics. This result holds because the sound-proof approximations considered in this paper only alter the form of $\Sigma$ but not its dependence on the core $(\eta, S, p, \Phi)$ variables. This result is important, because it demonstrates that it is possible to formulate sound-proof approximations in a way that makes their energetics and thermodynamics parallel that of the fully compressible Navier-Stokes equations, whose feasibility had remained unclear until now. For instance, Huang et al. (2001) argued that the Boussinesq equations are fundamentally affected by significant energy errors, but the present results establish that this can easily be avoided by adopting the above formalism and the proposed modernisation of the Boussinesq system.

This work underscores the fundamental importance of static energy for elucidating the energetics and thermodynamics of sound-proof approximations to the equations of motion. Physically, this is because the static energy encapsulates all fundamental information about the potential energy and pressure-geopotential gradient force that it generates, so that energetically and thermodynamically consistent approximations can be simply constructed by approximating the static energy to various orders of accuracy while respecting a few simple rules. From a practical standpoint, this approach naturally gives rise to two well-known important types of consistent approximations: pseudo-incompressible (PI) and (modernised) anelastic (AN), but which are here re-derived more transparently and succinctly while also clarifying some aspects that are not generally discussed, such as that pertaining to the choice of reference states.

The proposed approximation procedure is similar in many respects to Hamilton's Principle Asymptotics (HPA) advocated by Holm et al. (2002) and leveraged by many authors, e.g. Vasil et al. (2013), Tort and Dubos (2014b,a), Dubos and Voitus (2014), Dubos (2018). Holm et al. (2002) argue that approximations should be developed directly at the level of the Lagrangian from which the adiabatic equations of motion derive via Hamilton's least action principle, because this guarantees the existence of conservation laws whenever they should hold, i.e. when the Lagrangian presents the required symmetries. For the adiabatic part of the equations of motions, our method is thus a special case of HPA, limited to approximations of the Lagrangian that modify only static energy, thus motivating the name Static Energy Asymptotics (SEA) used in this paper to refer to our method. While more restrictive than the full HPA, SEA is also more comprehensive in that it deals consistently with diabatic processes. SEA is thus a stripped-down version of HPA applied to variational principles including diabatic processes, circumventing the advanced mathematical machinery used by Eldred and Gay-Balmaz (2021).

In the context of numerical ocean modelling, the PI and AN approximations are both superior to the standard SBA; their potential energy and thermodynamic potentials are nearly identical to their exact counterparts so that their energetics is more easily comparable to that of the fully compressible equations and do not require the introduction of artificial and ad-hoc thermodynamic potentials. The PI and AN approximations can therefore leverage the full power the Gibbs Sea Water (GSW) library developed as part of the new TEOS-10 equation of state (IOC, SCOR, IAPSO, 2010; Pawlowicz et al., 2012). Both approximations only affect the local APE density by removing its available
compressible part $Π_1$ but leave the local BPE unaffected. In contrast, how to define the local BPE for the SBA is much less clear, given that closing its energetics requires the construction of ad-hoc potentials. In relation to previous work, our seawater PI approximation (which is essentially a special case of the “atmospheric” PI approximation with seawater thermodynamics) provides a physically more transparent and more systematic foundation for the PI models that Dewar et al. (2016) proposed to improve on existing SBA-based OGCMs. As to our seawater AN approximation, it unifies the modernised-Boussinesq and anelastic approximations into a single framework, some of its elements being found in previous works by Pauluis (2008), Young (2010), Tailleux (2012). Importantly, it shows that to construct a seawater anelastic approximation, linearisation of the equation of state is not necessary, as previously established by Pauluis (2008) but in contrast to what was assumed in Ingersoll (2005). From a practical viewpoint, implementations of the PI models would require non-trivial alterations to existing OGCM codes and remains to be attempted. In contrast, the seawater generalised AN approximation has essentially the same mathematical structure as that of the standard SBA, so that its practical implementation should be straightforward. As doing so would significantly improve the energetics and thermodynamics of existing OGCMs at little to no additional cost, we recommend that efforts should be devoted in the future to adopt it. A natural choice of reference pressure and density profiles would be the analytical profiles entering the formulation of the analytical form of thermodynamic neutral density proposed in Tailleux (2021).

We believe that our findings can demystify a traditionally intricate and obscure aspect of atmospheric and ocean modelling. The issue of whether an approximation to the equations of motion is energetically and thermodynamically consistent has often been perceived as arcane, lacking systematic and general rules. We hope that our work will help atmospheric and ocean modellers realise that this issue is simpler than previously assumed, and that our findings will contribute to enhancing the physical realism of atmospheric and oceanic models, thus addressing the concerns raised by Lauritzen et al. (2022).

### CRediT authorship contribution statement

**Rémi Tailleux:** Writing – original draft, Methodology, Investigation, Formal analysis, Conceptualization. **Thomas Dubos:** Writing – review & editing, Methodology, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

No data was used for the research described in the article.

### Acknowledgements

This work was supported by an International Exchange award from the Royal Society number IES\R3\223161. This work is part of the AWACA project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 951596). The authors thank four anonymous reviewers for their useful and thorough remarks, which greatly helped improve the final manuscript.

### References

Alberty, R.A., 1994. Legendre transforms in chemical thermodynamics. Chem. Rev. 94, 1457-1482. http://dx.doi.org/10.1021/cr00030a001.

Cotter, C.J., Holm, D.D., 2014. Variational formulations of sound-proof models. J. R. Meteorol. Soc. 140, 1966-1973. http://dx.doi.org/10.1002/qj.2260.

Crocco, L., 1937. Eine neue Stromfunktion für die Erforschung der Bewegung der Gase mit Rotation. Z. Angew. Math. Mech. 17, 1-7. http://dx.doi.org/10.1002/zamm.19370170103.

de Szoeke, R.A., 2000. Equations of motion using thermodynamic coordinates. J. Phys. Oceanogr. 30, 2184-2829.

de Szoeke, R.A., Samelson, R.M., 2002. The duality between the Boussinesq and non-Boussinesq equations of motion. J. Phys. Oceanogr. 32, 2194-2203.

Dewar, W.K., Schoonover, J., McDougall, T.J., Klein, R., 2016. Semicompressible ocean thermodynamics and Boussinesq energy conservation. Fluids 1.

Dubos, T., 2018. Covariant structure of models of geophysical fluid motion. Phys. Rev. Lett. 120, 034501. http://dx.doi.org/10.1103/PhysRevLett.120.034501.

Dubos, T., Voitus, F., 2014. A semihydrostatic theory of gravity-dominated compressible flow. J. Atmos. Sci. 71, 4621-4638. http://dx.doi.org/10.1175/JAS-D-14-0080.1.

Durran, D.R., 1989. Improving the anelastic approximation. J. Atmos. Sci. 46, 1453-1461.

Durran, D.R., 2008. A physically motivated approach for filtering acoustic waves from the equations governing compressile stratified flows. J. Fluid Mech. 601, 365-379.

Eden, C., 2015. Revisiting the energetics of the ocean in Boussinesq approximation. J. Phys. Oceanogr. 45, 630-637.

Eden, C., 2016. Closing the energy cycle in an ocean model. Ocean Model. 101, 30-42.

Eden, C., Czeschel, L., Olbers, D., 2014. Toward energetically consistent ocean models. J. Phys. Oceanogr. 44, 3160-3184. http://dx.doi.org/10.1175/JPO-D-13-0260.1.

Eldred, C., Gay-Balmaz, F., 2021. Thermodynamically consistent semi-compressible fluids: a variational perspective. J. Phys. A 54, 345701. http://dx.doi.org/10.1088/1751-8121/ac1384.

Eldred, C., Taylor, M., Guba, O., 2022. Thermodynamically consistent versions of approximations used in modelling moist air. Q. J. R. Meteorol. Soc. 148, 3184-3210.

Holm, D.D., Marsden, J.E., Ratiu, T.S., 2002. The Euler-Poincaré equations in geophysical fluid dynamics. In: Norbury, I. (Ed.), Large-Scale Atmosphere-Ocean Dynamics 2: Geometric Methods and Models. Cambridge University Press, Cambridge, pp. 251-299.

Huang, R.X., Xiangze, J., Xuehong, Z., 2001. An oceanic general circulation model in pressure coordinates. Adv. Atm. Sci. 18, 1-22.

Ingersoll, A.P., 2005. Boussinesq and Anelastic approximations revisited: potential energy release during thermobaric instability. J. Phys. Oceanogr. 35, 1359-1369.

IOC, SCOR, IAPSO, 2010. The international thermodynamic equation of seawater - 2010: Calculation and use of thermodynamic properties. In: Manual and Guides No. 56. Intergovernmental Oceanographic Commission, UNESCO (English).

Klein, R., Pauluis, O., 2012. Thermodynamic consistency of a pseudoincompressible approximation for general equations of state. J. Atmos. Sci. 69, 961-968.

Lauritzen, P.H., Kevlahan, N.K.R., Toniazzo, T., Eldred, C., Dubos, T., Gassmann, A., Larson, V.E., Jablonowski, C., Guba, O., Shipway, B., Harrop, B.E., Lemarié, F., Tailleux, R., Herrington, A.R., Large, W., Rasch, P.J., Donahue, A.S., Wan, H., Conley, A., Bacmeister, J.T., 2022. Reconciling and improving formulations for thermodynamics and conservation principles in Earth System Models (ESMs). J. Adv. Model. Earth Syst. 14, http://dx.doi.org/10.1029/2022MS003117.

Lipps, F.B., Hemler, R.S., 1982. A scale analysis of deep moist convection and some related numerical calculations. J. Atmos. Sci. 39, 2192-2210.

Lorenz, E.N., 1955. Available potential energy and the maintenance of the general circulation. Tellus 7, 157-167.

Losch, M., Adcroft, A., Campin, J.M., 2004. How sensitive are coarse general circulation models to fundamental approximations in the equations of motion? J. Phys. Oceanogr. 34, 306-3119.

McDougall, T.J., 1987. Thermobaricity, cabbeling, and water mass-conversion. J. Geophys. Res. Oceans 92, 5448-5464.

McDougall, T.J., 2003. Potential enthalpy: a conservative oceanic variable for evaluating heat content and heat fluxes. J. Phys. Oceanogr. 33, 945-963.

McDougall, T.J., Greatbatch, R.J., Lu, Y., 2002. On conservation equations in oceanography: how accurate are Boussinesq ocean models? J. Phys. Oceanogr. 32, 1574-1584.

Nycander, J., 2010. Horizontal convection with a non-linear equation of state: generalization of a theorem of Paparella and Young. Tellus 62A, 134-237.

Nycander, J., Hieronymus, M., Roquet, F., 2015. The nonlinear equation of state of sea water and the global water mass distribution. Geophys. Res. Lett. 42, 7714-7721. http://dx.doi.org/10.1002/2015GL065525.

Ogura, Y., Phillips, N.A., 1962. Scale analysis of deep and shallow convection in the atmosphere. J. Atmos. Sci. 19, 173-179.

Pauluis, O., 2008. Thermodynamic consistency of the Anelastic approximation for a moist atmosphere. J. Atmos. Sci. 65, 2719-2729.

Pawlowicz, R., McDougall, T.J., Feistel, R., Tailleux, R., 2012. An historical perspective on the development of the thermodynamic equation of seawater - 2010. Ocean Sci. 8, 161-174. http://dx.doi.org/10.5194/os-8-161-2012.
Saenz, J.A., Tailleux, R., Butler, E.D., Hughes, G.O., Oliver, K.I.C., 2015. Estimating lorenz's reference state in an ocean with a nonlinear equation of state for seawater. *J. Phys. Oceanogr.* 45, 1242–1257. http://dx.doi.org/10.1175/JPO-D-14-0105.1

Scinocca, J.F., Shepherd, T.G., 1992. Nonlinear wave-activity conservation laws and hamiltonian structure for the two-dimensional Anelastic equations. *J. Atmos. Sci.* 49, 5–28.

Shepherd, T.G., 1993. A unified theory of available potential energy. *Atmos.-Ocean* 31, 1–26.

Tailleux, R., 2010a. Entropy versus APE production: on the buoyancy power input in the oceans energy cycle. *Geophys. Res. Lett.* 37, http://dx.doi.org/10.1029/2010GL044962.

Tailleux, R., 2010b. Identifying and quantifying nonconservative energy production/destruction in hydrostatic Boussinesq primitive equation models. *Ocean Model.* 34, 125–136. http://dx.doi.org/10.1016/j.ocemod.2010.05.003.

Tailleux, R., 2012. Thermodynamics/dynamics coupling in weakly compressible turbulent stratified fluids. *ISRN Thermodyn.* 2012, 609701. http://dx.doi.org/10.5402/2012/609701.

Tailleux, R., 2013. Available potential energy density for a multicomponent Boussinesq fluid with a nonlinear equation of state. *J. Fluid Mech.* 735, 499–518. http://dx.doi.org/10.1017/jfm.2013.509.

Tailleux, R., 2018. Local available energetics of multicomponent compressible stratified fluids. *J. Fluid Mech.* 842, http://dx.doi.org/10.1017/jfm.2018.196.

Tailleux, R., 2021. Spiciness theory revisited, with new views on neutral density, orthogonality and passiveness. *Ocean Sci.* 17, 203–219. http://dx.doi.org/10.5194/os-17-203-2021.

Tort, M., Dubos, T., 2014a. Dynamically consistent shallow-atmosphere equations with a complete Coriolis force. *Q. J. R. Meteorol. Soc.* 140, 2388–2392. http://dx.doi.org/10.1002/qj.2274.

Tort, M., Dubos, T., 2014b. Usual approximations to the equations of atmospheric motion: A variational perspective. *J. Atmos. Sci.* 71, 2452–2466. http://dx.doi.org/10.1175/JAS-D-13-0339.1.

Vasil, G.M., Lecoanet, D., Brown, B.P., Wood, T.S., Zweibel, E.G., 2013. Energy conservation and gravity waves in sound-proof treatment of stellar interiors. *Part II. Lagrangian constrained analysis. Astrophys. J.* 773, 1–23.

Vazsonyi, A., 1945. On rotational gas flows. *Q. Appl. Math.* 3, 29–37.

Woods, L.C., 1975. The Thermodynamics of Fluid Systems. Clarendon Press.

Young, W.R., 2010. Dynamic enthalpy, conservative temperature, and the seawater boussinesq approximation. *J. Phys. Oceanogr.* 40, 394–400.

Young, W.R., 2012. An exact thickness-weighted average formulation of the Boussinesq equations. *J. Phys. Oceanogr.* 692–707.