A Simple and Transparent Method for Improving the Energetics and Thermodynamics of Seawater Approximations: Static Energy Asymptotics (SEA)

The static energy encodes all possible information about the thermodynamics and potential energy (and all related forces) of stratiﬁed geophysical ﬂuids. 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 a the anelastic (AN) approximation. In all cases, the background and available parts of the potential energy (in Lorenz sense) can be constructed to coincide exactly and approximately with their exact counterparts respectively, and in both cases be expressed in terms of the exact (as opposed to ad-hoc) thermodynamic potentials recently made available by TEOS-10. 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 numerical ocean models. 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 TEOS-10 Gibbs Sea Water (GSW) library.


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 et al. (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 soundproofing 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 v ρ = ρ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 ∇ • v ρ = 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 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 subgridscale energy reservoirs, e.g., Eden et al. (2014); Eden (2015Eden ( , 2016)).
Following relatively recent progress (Young, 2010;Tailleux, 2010bTailleux, , 2012;;Eden et al., 2014;Eden, 2015Eden, , 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 et al., 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 accomodate 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 (1989Durran ( , 2008a) ) 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 formulation 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 (2015Eden ( , 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 sci-entists, 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 and background potential energy, as per Tailleux (2018) local formulation of Lorenz (1955) theory of available potential energy.Section 7 summarises and discusses the results.

Energy conservation and structure of static energy
To understand how to construct energetically and thermodynamically consistent, 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.

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 Here, h = h(η, S, p) is the specific enthalpy, v = (u, v, w) is the threedimensional velocity field, ρ is density, p is pressure, Φ = gz is the geopotential, g is the acceleration of gravity, z is height, S is salinity, η is specific entropy, with J S and J η the molecular fluxes of salt and entropy respectively.The term ηirr > 0 represents the irreversible rate of entropy production, which the second law of thermodynamics requires to be positive definite.Finally, F = ρ −1 ∇•S is the viscous force, which may be related to the deviatoric stress tensor S.
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 η and J S , as well as with the constraint that ηirr > 0 be positive definite.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 = contant and relative chemical potential µ = constant.The simplest way to achieve this is via assuming that the molecular fluxes of η and S be linearly related to the generalised forces ∇T /T and ∇µ/T as follows with L ηη , L ss , L ηs and L sη to be determined empirically from laboratory measurements.Energy conservation can be shown to impose that the following quantity be equal to the divergence of an appropriately define energy flux J E .This in turn constrains the non-viscous irreversible entropy production to take the form which is positive definite provided that Finally, the Onsager reciprocity relationship allows for the simplification L ηs = L sη .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 ( 4) 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 build to satisfy the law of energy conservation, the conserved total energy being 2 the kinetic energy (KE), and the potential energy (PE), itself the sum of gravitational potential energy Φ(z) = gz and internal energy e = h − p/ρ.Note that the Navier-Stokes equations contain only information about the three partial derivatives of h, meaning that they are only sufficient to recover h up to an irrelevant integration constant.(Note though that while ν and T are unambiguously defined, µ 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 et al. ( 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 ( η = Ṡ = 0), for which knowledge of the thermohaline derivatives of h (T and µ) is no longer formally needed for integrating the equations forward in time.In this case, the Navier-Stokes equations reduce to the Euler equations and only contain enough information to recover h (or e) up to an indeterminate function of (η, 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 (η, 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 θ or Conservative Temperature Θ in place of entropy as h = ĥ(S, θ, p) or h(S, Θ, p), as is usually preferred in oceanographic practice.To avoid losing information in that case, it becomes necessary to also supply the passage relations η = η(S, θ, p) or η = η(S, Θ, p) allowing one to reconstruct h in terms of its canonical variables, e.g., IOC et al. (2010).This is needed, for instance, to compute variables such as in-situ temperature T , whose expression is which clearly requires knowledge of both ĥ and η.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 recovery strategies.

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 where ρ b is the constant Boussinesq reference density, which in the NEMO OGCM is generally chosen to be Although (15-20) has formed the basis for a majority of OGCMs, that they admit an energy conservation principle was only clarified very recently by Young (2010), Tailleux (2012), andEden (2015) among others, building upon previous ideas by Ingersoll (2005) and Pauluis (2008).Starting from the kinetic energy equation, obtained by taking the scalar product of ( 15) with u combined with ( 16) multipled 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 where the Boussinesq pseudo potential energy E p,bou and the thermodynamic efficiencies Υ S and Υ θ are defined by where z ⋆ is an arbitrary reference depth.Eq. ( 22) defines an energy conservation principle for u 2 /2 + E p,bou for inviscid motions in the absence of diabatic sources/sinks of θ 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.

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 ρ −1 ∇p+∇Φ, potential energy E p , and dynamics/thermodynamics coupling in terms of the static energy by means of the following relations: where in the case of a fully compressible fluid, the partial derivatives of 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 may also be written in thermodynamic form as follows: 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 Σ F C 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) Σ F C is some general function of (η, S, p, Φ); 2) that the density entering the continuity equation ( 2) relates to the pressure derivative of Σ F C via Thus, taking the scalar product of (32) with v, using the facts that by definition v • ∇C = DC/Dt − ∂C/∂t = Ċ − ∂C/∂t for any scalar field C(x, t) and that ∂Φ/∂t = 0, yields Next, multplying the result by (∂Σ F C /∂p) −1 = ρ, accounting for the continuity equation (2), yields Eq. ( 35) makes it clear that E total is conserved in the absence of viscous (F = 0) and diabatic ( η = Ṡ = 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., Physically, Eq. ( 37) in turn constrains the form of ηirr in terms of the molecular fluxes J η and J s as well as with the gradients of ∇T and ∇µ 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 where is the standard Bernoulli function attached to Σ F C .Eq. ( 38) is the appropriate form of the energy conservation principle for establishing the Bernoulli theorem, namely that B F C is conserved along adiabatic, isohaline, inviscid, steady fluid parcel trajectories.

Pseudo-compressible approximation
To summarise, we established that under the conditions stated above, the equations of motion formulated in terms of the static energy Σ F C (η, S, p, Φ) as per (31) or (32), together with the continuity equation ( 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 regardless of the exact form of Σ F C .As a consequence, the same result must hold for any approximation to the equations of motion obtained by any structure-preserving approximation of Σ F C , 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 pseudoincompressible (PI) approximation (Durran, 1989(Durran, , 2008b;;Dewar et al., 2016;Eldred et al., 2022) via replacing Σ F C by an approximation Σ P I defined in terms of the leading order terms of a Taylor series expansion around the reference pressure p R (Φ) as follows: where δp = p − p R (Φ) is the assumed small pressure anomaly and R P I = O(δp 2 ) the residual term.Importantly, Eq. ( 41) shows that Σ P I = Σ P I (η, S, p, Φ) retains the same functional form as Σ F C and hence that the above results must apply.Evaluating the (η, S, p, Φ) derivatives of Σ P I yields and show that these remain close to the derivatives of Σ F C , with relative minor impact on ρ, T , and µ, the main effect being the actual pressure p being replaced by the reference pressure p R (Φ) and the presence of a small correction O(δp), where ρ R (Φ) = −dp R /dΦ denotes the reference density in hydrostatic equilibrium with p R .As a result, the above approximation affects the Crocco theorem, momentum balance equation, and continuity equation as follows: noting that, per (33) applied to Σ P I instead of Σ F C , ρ must also be replaced by ρ ⋆ in the tracers equations for specific entropy and salinity (3) and (4).
Since Σ P I has the same functional structure as Σ F C , 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,P I , with provided that the diabatic terms η and Ṡ are constrained to satisfy In that case, (49) shows that the approximation to the potential energy E p,P I 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 , µ, and ρ by T P I , µ P I , and ρ ⋆ in the expressions for the diabatic terms η and Ṡ, as previously found by Klein and Pauluis (2012) and Eldred and Gay-Balmaz (2021).

Remarks on the choice of p R (Φ)
Physically, the accuracy of the PI approximation is determined by the magnitude of the residual R P I in (41), which may be rewritten in the following equivalent forms: Because υ p = −1/(ρ 2 c 2 s ) < 0 where c s is the speed of sound, it follows that R P I is also negative.Moreover, it is also easily shown that R P I is quadratic in δp at leading order, viz., and hence that its magnitude is directly controlled by the choice of p R (Φ) and its distance to p.The reader familiar with the local theory APE will may recognise that R 1 is approximately equal to the opposite of available compressible energy (ACE) denoted by Π 1 in Tailleux (2018).As a consequence, which establishes, perhaps counter-intuitively, that Σ P I represents an overestimate of the true Σ F C .In the literature, how to specify p R (Φ) is rarely is every discussed.The fact that R P I is sign definite suggests, however, that p R (Φ) should be defined to minimise the volume integral of R P I .In the case where the denominator ρ 2 ⋆ c 2 s⋆ can be treated as approximately constant, this would amount to define p R (Φ) as the horizontal-mean pressure.

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 it is important to understand and be aware of for correctly interpreting some results of the literature.
As shown below, a alternative energy conservation principle may be shown to exist in the PI approximation because (from Eq. ( 41)) can alternatively be interpreted as a function of (η, S, δp, Φ) (denoted with a hat).Since δp = p − p R (Φ) is a function of p and Φ only, the pressure-geopotential gradient force can therefore be written equivalently as the (δp, Φ) partial derivatives of ΣP I being given by while in contrast its thermodynamic derivatives remain unaffected From ( 55), it follows that the PI momentum balance equations may alternatively be written in the form which from the same structural arguments developed previously implies the existence of a degenerate form of energy conservation for the pseudo total energy E total = E k + E p,pseudo , with so that E p,pseudo and E p,P I differ, specifically by the term p R (z)/ρ ⋆ .However, as the latter quantity may be shown to satisfy the conservation law it follows that the conservation principle for E k +E p,pseudo simply results from adding (62) to the conservation principle for the correct energy E k + E p,P I .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 fullycompressible energy, a point that many authors do not appear to be aware of, e.g., see Eq. ( 26) of Dewar et al. (2016) for instance.
The AN approximation Σ AN retains the same (η, S, p, Φ) functional dependence characterising Σ F C and Σ P I .Its partial derivatives are: (with υ ‡ = ρ −1 ‡ ) and can be seen to remain close to the derivatives of Σ P I , albeit with a loss of accuracy due to the absence of the small correction O(δp) and to ρ ⋆ = ρ ‡ .
The fact that Σ AN has the same structure as Σ F C or Σ P I means that provided that its diabatic and viscous terms satisfy the constraint analogous to (9), it must conserve the total energy E k + E p,AN , with the AN potential energy being given by (69) Eq. ( 69) shows that E p,AN differs only from E p,P I by the the energy error term ∆E p,spurious = p R (υ ⋆ − υ ‡ ).This error can be minimised by defining υ ‡ as the temporally and horizontally averaged value of υ ⋆ for instance.
Eqs. (64)(65)(66)(67) show that the CVT can be written in the standard form However, a more convenient and concise expression of the CVT can be obtained by taking advantage of the possibility of regarding Σ AN = ΣAN (η, S, δp/ρ ‡ , Φ) as a function of (η, S, δp/ρ ‡ , Φ) instead, for which the (η, S, Φ) derivatives can be shown to be where b AN is recognised as the 'buoyancy' defined relative to ρ R (Φ), which differs from the traditional Boussinesq buoyancy as further discussed below.It follows that the CVT (70) may be rewritten in the form Eq. ( 74) differs from ( 70) in that the density ρ ‡ is now absorbed within the anomalous pressure gradient term, while only the buoyancy b AN multiplies ∇Φ, thus making it possible to rewrite the standard form of momentum balance in the form 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 ρ ‡ depends on z, as further discussed in the next section.

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 Φ, which is more standard in oceanographic practice) where ρ ‡ (z) does not need to coincide with ρ 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) = −ρ b gz, with ρ b = constant, arguments developed in next section suggest that p R (z) and ρ 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 the legacy SBA model (15)(16)(17)(18)(19)(20), and hence that it can be solved algorithmically and numerically in exactly the same way regardless of whether ρ ‡ 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 nonapproximated 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., 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 ρ R (z) that it uses is much clearer than that for specifying the constant Boussinesq reference density ρ 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 ρ ‡ could also be chosen as a constant reference Boussinesq density ρ b , this simplification offers no advantage over using a depth-dependent ρ 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 modernized version proposed here and to the anelastic system.

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, as justified by Lorenz (1955) theory of available potential energy (APE) for instance.Since these two flavours behave dynamically very differently (Tailleux, 2010b), it follows that any discussion of energetics and thermodynamics consistency cannot be complete without understanding how sound-proof approximations individually affect the available and non-available parts of the 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 (Σ dyn ) and passive (Σ heat ) components as follows: Since Σ heat is independent of p and Φ, it follows that only Σ dyn is dynamically relevant for predicting the pressure-geopotential gradient force entering the standard form of momentum balance equations, thus justifying the dynamically irrelevant character of Σ heat .In the following, we briefly discuss the two main approaches for defining Σ dyn and Σ 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 Σ dyn should only require knowledge of ∂Σ/∂p and ∂Σ/∂Φ (regardless of how Σ 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.

Available potential energy (APE) theory
From a theoretical viewpoint, the most rigorous and physically-based partitioning of Σ F C is arguably the one rooted in APE theory.In this approach, Σ dyn is defined as the departure of Σ F C 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 ρ 0 (z), that is, where z R = z R (η, S) is the fluid parcel's reference position.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(Tailleux, , 2018) ) ρ(η, S, p 0 (z R )) = ρ 0 (z R ). ( 85) With Σ heat defined as per (84), we have where Π 1 and Π 2 are so-called available compressible energy (ACE) density and APE density respectively, Π 2 = h(η, S, p 0 (z)) − h(η, S, p 0 (z R )) + g(z − z R ). (88) As shown by Tailleux (2018), these may be written as where b(η, S, z) = −g[1 − ρ 0 (z)υ(η, S, p 0 (z))] is the buoyancy defined relative to the reference density ρ 0 (z).This establishes, therefore, that Π 1 , Π 2 , and therefore Σ dyn can be constructed only from the knowledge of υ(η, S, p), g, as well as from one state of the system from which p 0 (z) and ρ 0 (z) are constructed.
Let us now seek to understand how the AN approximation affects the dynamics/thermodynamics partitioning of Σ AN , By considering the value of Σ AN in the Lorenz background reference state, in which z = z R and p = p 0 (z R ), it is easily realised that Σ 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 Σ AN becomes where the expression for Π 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 Π 1 , while also approximating the compressible work term (p − p 0 (z))/ρ by (p − p 0 (z))/ρ ‡ .To avoid the proliferation of reference density profiles, the most logical choice is to use ρ ‡ = ρ 0 (z), which provides a natural solution of how to choose ρ ‡ .

Potential/Dynamic enthalpy partitioning
We briefly remark that another kind of partitioning is based on defining Σ heat as the 'potential' value of Σ F C referenced to the ocean surface at z = 0, viz., Σ heat = h(η, S, p R (0)) + Φ(0), ( 93) (de Szoeke, 2000;Young, 2010;Nycander, 2010) where p R (z) is an arbitrary hydrostatic reference pressure field defining the reference density field via ρ R (z) = −p ′ R (z)/g.If p R (0) = p a and Φ(0) = 0, where p a is the mean surface atmospheric pressure, Σ heat then coincides with McDougall (2003)'s potential enthalpy, the currently recommended variable for defining heat content in the oceans following the adoption of TEOS-10 (IOC et al., 2010;Pawlowicz et al., 2012).With Σ heat defined as per (93), the expression for Σ dyn can easily be verified to be with Π 1 = h(η, S, p) − h(η, S, p R (z)) + (p R (z) − p)/p ≥ 0 the ACE defined in terms of p R (z) rather than p 0 (z), and h ‡ the so-called dynamic enthalpy, defined by ] the buoyancy defined relative to the reference density ρ R (z).It is easily seen that ( 94) and ( 93) 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 Σ AN , so that Σ AN,heat = Σ heat .The AN approximation only affects the dynamical component by getting rid of Π 1 in (94), so that The potential/dynamic enthalpy based partitioning of Σ AN , however, lacks any of the advantageous attributes and clear physical interpretation of the BPE/APE based partitioning; its usefulness is therefore questionable.

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 soundproof approximations of such equations, under the generic form Σ = Σ(η, S, p, Φ) (97) T where E tot = v 2 /2 + Σ − p∂Σ/∂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 Σ but not its dependence on the core (η, S, p, Φ) variables.This result is important, because it demonstrates that it is possible to formulate sound-proof approximations in a way that make 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 modernization 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 Halmilton'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 at 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 fully compressible 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 et al., 2010;Pawlowicz et al., 2012).Both approximations only affect the local APE density by removing its available compressible part Π 1 but leaves the local BPE unaffected.In contrast, how to define the local BPE for the SBA is much less clear, giving 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 modernized-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 is 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).
search Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 951596)

TS
J s = −L ηs ∇T − L ss ∇µ (105) As shown in this paper, such equations conserve total mass defined as and total energy defined as the volume integral of ∂Σ ∂p −1 E tot ,