# Energetically consistent localised APE budgets for local and regional studies of stratified flow energetics

Rémi Tailleux, Guillaume Roullet

► To cite this version:

Rémi Tailleux, Guillaume Roullet. Energetically consistent localised APE budgets for local and regional studies of stratified flow energetics. Ocean Modelling, 2025, 197 (2), pp.102579. 10.1016/j.ocemod.2025.102579. hal-05322754

HAL Id: hal-05322754

https://hal.univ-brest.fr/hal-05322754v1

Submitted on 20 Oct 2025

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
# Energetically consistent localised APE budgets for local and regional studies of stratified flow energetics

Rémi Tailoux<sup>a, ID*</sup>, Guillaume Roullet<sup>b</sup>

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

## ARTICLE INFO

### Keywords:

- Available potential energy

- Energy cascades

- Energetically consistent ocean models

- Turbulent mixing parameterisations

## ABSTRACT

Because it allows a rigorous separation between reversible and irreversible processes, the concept of available potential energy (APE) has become central to the study of turbulent stratified fluids. In ocean modelling, it is fundamental to the parameterisation of meso-scale ocean eddies and of the turbulent mixing of heat and salt. However, how to apply APE theory consistently to local or regional subdomains has been a longstanding source of confusion due to the globally defined Lorenz reference state entering the definition of APE and of buoyancy forces being generally thought to be meaningless in those cases. In practice, this is often remedied by introducing heuristic 'localised' forms of APE density depending uniquely on region-specific reference states, possibly diverging significantly from the global Lorenz reference state. In this paper, we argue that across-scale energy transfers can only be consistently described if localised forms of APE density are defined as the eddy APE component of an exact mean/eddy decomposition of the APE density, for which a new physically more intuitive and mathematically simpler framework is proposed. The eddy APE density thus defined exhibits a much weaker dependency on the global Lorenz reference state than the mean APE, in agreement with physical intuition, but with a different structure than that of existing heuristic localised APE forms. Our framework establishes a rigorous physical basis for linking parameterised energy transfers to molecular viscous and diffusive dissipation rates. We illustrate its potential usefulness by discussing the energetics implications of standard advective and diffusive parameterisations of the turbulent density flux, which reveals potential new sources of numerical instability in ocean models.

## 1. Introduction

Turbulent stratified flows exhibit a complex interplay between reversible and irreversible processes, each providing crucial insights into the other. Reversible aspects are typically associated with adiabatic stirring, which deforms isopycnal surfaces, increases their areas, and magnifies irreversible effects by enhancing tracer gradients. This process leads to the dissipation of mechanical energy and tracer variances at increasingly smaller scales through molecular and diffusive processes (Eckart, 1948). The concept of available potential energy (APE), originally formulated by Margules (1903) and Lorenz (1955) and later adapted to the study of turbulent stratified mixing by Winters et al. (1995), serves as a key tool for distinguishing between reversible and irreversible processes. APE theory posits that the potential energy (PE) of any stratified fluid can be partitioned into a component (the APE) available for reversible conversions with kinetic energy (KE), and a component (the background potential energy, BPE) that is not. In

Lorenz's approach, the BPE is defined as the PE of a flattened state of minimum potential energy obtainable from the actual state through an adiabatic rearrangement of mass (Lorenz, 1955). Consequently, APE theory provides a natural framework for assigning distinct energetic signatures to reversible and irreversible processes (Butler et al., 2013). Reversible processes affect the APE of the fluid while leaving the BPE unaffected, whereas irreversible processes entail an energy transfer between the APE and BPE. In most cases, the net transfer occurs from the APE to the BPE; however, the reverse conversion is occasionally possible, as observed in double-diffusive instabilities (e.g., Middleton and Taylor, 2020; Middleton et al., 2021; Tailoux, 2024). In the local theory of APE, the conversion rate between APE and BPE is generally referred to as the APE dissipation rate, denoted as $\epsilon_p$. Although initially introduced in the context of Boussinesq fluids, the nature of $\epsilon_p$ in the fully compressible Navier-Stokes equations has been discussed

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

<footnote><sup>b</sup> Université de Bretagne Occidentale, CNRS, IRD, Ifremer, Laboratoire d'Océanographie Physique et Spatiale (LOPS), rue Dumont Durville, Plouzané, 29280, France</footnote>

<footnote>* Corresponding author.
E-mail address: r.g.j.tailleux@reading.ac.uk (R. Tailoux).</footnote>
by Tailleux (2009, 2013c) and Tailleux (2024). In the global APE framework of Winters et al. (1995), the volume-integrated APE dissipation rate is represented by the term $\Phi_d - \Phi_l$. Together with the viscous dissipation rate $\epsilon_k$, the sum $\epsilon_k + \epsilon_p$ represents the total dissipation of mechanical energy, defined as the sum of KE and APE.

While recent discussions of APE have primarily been framed within the global APE framework of Winters et al. (1995), the concept's importance and utility were first recognised by Oakey (1982) and Gargett and Holloway (1984) through the derivation of a local APE budget for a quadratic, non-negative APE density. This was achieved by rescaling the budget of temperature variance and linking it to the mechanical energy budget. This result was significant for connecting turbulent mixing to the mechanical energy budget, as it introduced a conversion term with kinetic energy and related the APE dissipation rate $\epsilon_p$ to the dissipation of temperature variance:

$$ \epsilon_p = \frac{g\alpha\kappa_T |\nabla\theta'|^2}{d\bar{\theta}/dz}. \qquad (1) $$

In (1), g represents the acceleration due to gravity, a is the thermal expansion coefficient, $\kappa_T$ is the thermal diffusivity, and $\theta$ is the potential temperature with mean $\bar{\theta}$ and perturbation $\theta'$. The relative importance of diffusive and viscous effects in dissipating mechanical energy can be quantified using the dissipation ratio $\Gamma = \epsilon_p/\epsilon_k$, a commonly used measure of mixing efficiency, often considered to be close to 0.2. Both $\epsilon_p$ and $\epsilon_k$ can be used to define the turbulent diapycnal mixing according to the formula

$$ K_p = \frac{\epsilon_p}{\bar{N}^2} = \frac{\Gamma \epsilon_k}{\bar{N}^2}, \qquad (2) $$

(e.g., Lindborg, 2006). Over the past three decades, the somewhat ad-hoc approach to APE developed by Oakey (1982) and Gargett and Holloway (1984) has been superseded by the exact global APE framework of Winters et al. (1995) and the exact finite-amplitude local APE framework first developed by Andrews (1981) and Holliday and McIntyre (1981), and subsequently extended by Shepherd (1993), Scotti et al. (2006), Roullet and Klein (2009), Tailleux (2013b), Scotti and White (2014), Zemskova et al. (2015), and Tailleux (2018), among others.

Reversible and irreversible effects can also be described in terms of energy transfers between different scales of motion. In the atmosphere, there has been extensive discussion about the physical explanation for the turbulent energy cascade affecting both kinetic and potential energy. Lindborg (2006) developed a theory suggesting a forward energy cascade with energy spectra given by

$$ E_{K_h} = C_1 \epsilon_K^{2/3} k_h^{-5/3}, \quad E_{p_h} = C_2 \epsilon_p k_h^{-5/3} \epsilon_K^{-1/3} \qquad (3) $$

with $C_1 \approx C_2$. Interestingly, this theory predicts that the ratio of potential energy to kinetic energy spectra at all scales is

$$ \frac{E_{p_h}}{E_{K_h}} \approx \frac{\epsilon_p}{\epsilon_K} = \Gamma, \qquad (4) $$

thus suggesting that the constraint imposed by the dissipation ratio $\Gamma$ extends to significantly larger scales than the dissipation scales. Intriguingly, the ratio $G_A/G_K$ of the APE production rate $G_A$ by surface buoyancy fluxes to the wind power input into geostrophic motions $G_K$ also appears to be close to the dissipation ratio $\Gamma$. Indeed, based on Zemskova et al. (2015), $G_A \approx 0.5$ TW, and $G_K \approx 2-3$ TW, yielding a ratio $G_A/G_K \approx 0.17-0.25$. Strong constraints on the relative importance of the KE and APE energy spectra and energy transfers therefore appear to exist over a wide range of scales. Understanding the nature of such constraints and exploiting them to constrain the parameterisations of subgridscale processes controlling the unresolved KE and APE energy cascades forms the basis for the development of energetically consistent numerical ocean models.

Lorenz introduced the global APE framework to explain how the atmospheric circulation is maintained against dissipation. He introduced

the Lorenz Energy Cycle (LEC) as a means to achieve this, partitioning the KE and APE reservoirs into mean and eddy components. However, because the integrand defining the APE is not positive definite, Lorenz had to rely on various manipulations and integrations by parts to rewrite the volume-integrated APE as the volume integral of a positive definite quadratic quantity that could then be split into mean and eddy components. Setting aside compressible effects, the local APE framework defines the APE density of a fluid parcel in terms of the work against buoyancy forces relative to the reference density profile $\rho_0(z)$, characterising the globally defined Lorenz reference state of minimum potential energy. However, this reference state is rarely considered relevant for understanding the energetics of much smaller subdomains. Rather, it is generally implicitly assumed that buoyancy forces should be defined relative to some locally defined averaged density field characteristic of the local environment. To continue using APE theory in such cases, most existing approaches appear to be based on some heuristic 'localisation' of the local or global APE frameworks. A typical example of such localisation is the definition of eddy APE by Roullet et al. (2014), which physically modifies the exact form of local APE density

$$ E_a^{\text{exact}} = \frac{g}{\rho_*} \int_{z_0(\rho)}^{z} [\rho - \rho_0(\bar{z})] \, d\bar{z} \qquad (5) $$

into the following 'localised' form:

$$ E_a^{R14} = \frac{g}{\rho_*} \int_{z_m}^{z} [\rho - \rho_m(x, y, \bar{z})] \, d\bar{z}, \qquad (6) $$

with $\rho_m(x, y, z)$ representing some locally defined mean density field, and $z_m$ the level of neutral buoyancy satisfying $\rho_m(x, y, z_m) = \rho$, while $z_0(\rho)$ is such that $\rho_0(z_0(\rho)) = \rho$. While plausible and physically appealing, (6) is the source of much confusion in the literature about turbulent stratified mixing, as it no longer provides clarity on how to compute the reference state associated with the exact formula (5). Several studies have discussed the issue (e.g., Arthur et al., 2017; Wykes et al., 2015; Dewar and McWilliams, 2019) and found that different choices of reference states often lead to significantly different conclusions about the properties of mixing. To help visualise the nature of the problem, the left panel of Fig. 1 illustrates the differences between the exact (black) and heuristic (orange) APE densities predicted by (5) and (6) in a particular example for which $\rho_m$ and $\rho_0$ differ from each other, as is normally the case. The right panel anticipates the results of this paper, comparing the exact form of instantaneous eddy APE density resulting from an exact mean/eddy decomposition (black) with the heuristic APE density (orange, same as in left panel), associated with fluctuations around the mean density $\rho_m$. This schematic shows that while both approaches are associated with the same density anomalies $\rho'$, they are associated with different displacement anomalies $\zeta'$, potentially leading to significant differences.

Part of the difficulty or confusion surrounding this issue seems to arise from the insistence on discussing the energetics of individual mixing events independently of the energetics of the global ocean in which they are embedded. However, in a statistically steady state, the total KE+APE dissipation must balance the work done by the surface wind stress plus the APE production by the surface buoyancy fluxes. Considering that the APE production by surface buoyancy fluxes is always based on the exact APE density using the Lorenz reference density profile $\rho_0(z)$, it seems evident that if the APE dissipation rate based on (6) is sensitive to the choice of $\rho_m$, then there must be only one consistent way to define $\rho_m$ that can achieve the desired balance. The question addressed in this paper is: which one is it?

To ensure consistency between APE production at large scales and APE dissipation at molecular scales, the solution discussed in this paper is to define localised eddy APE densities in terms of the eddy component of an exact mean/eddy decomposition of the local APE density. As shown in this paper, this is especially important and necessary where the actual state departs significantly from the Lorenz reference state
<figure><img src="image_3.png" /><figcaption>Fig. 1. (a) Schematic illustrating the differences between the exact APE density (black) versus heuristic APE density (orange) as area under the curves, for a parcel with density $\rho$ at depth $z$; (b) Same as for left panel but for the instantaneous exact eddy APE density (black) versus heuristic APE density (orange) associated with density anomalies centred around the mean density $\rho_m(z)$ at depth $z$. Note that the black and orange areas can both be approximated as $\rho'\zeta'/2$ in terms of the same density anomaly $\rho'$ but different displacement anomaly $\zeta'$.</figcaption></figure>

and hence where the standard QG approximation becomes inaccurate. So far, a general theory for how to do so has been lacking, as existing results have been limited to a Boussinesq fluid with a linear equation of state (Scotti and White, 2014) or a dry atmosphere (Novak and Tailleux, 2018), using methods too specialised to be easily generalisable. It follows that to apply existing theory to a realistic ocean (e.g., Zemskova et al. (2015) or MacCready and Giddings (2016)), several approximations need to be made to circumvent the difficulties arising from the two-component nature of the equation of state and its thermobaric nonlinearity, thus diminishing its benefits. Perhaps because of this, most studies of the Lorenz energy cycle – while often acknowledging the existence of Scotti and White (2014) – still continue to rely on the standard QG approximation used by von Storch et al. (2012) and others.

In this paper, we revisit Scotti and White (2014)'s approach by making it more physically intuitive, mathematically simpler, and more easily generalisable to more complex fluids. To that end, we adopt as our starting point the local budget of APE density developed by Tailleux (2024), which more naturally generalises to two-component fluids. To summarise, our approach defines the mean APE density $E_a(\bar{\rho}, z)$ as the APE density of the mean density field $\bar{\rho}$, the eddy APE density as the residual $E_a^t = \bar{E}_a - E_a(\bar{\rho}, z)$, and proves the positive definite character of $E_a^t$ by linking it to the convexity of $E_a(\rho, z)$ with respect to $\rho$. The generality of such results suggests that for a realistic thermobaric ocean, the mean APE density should be defined as $E_a^m = E_a(\bar{S}, \bar{\theta}, z)$ and hence the eddy APE density as the residual $E_a^t = \bar{E}_a - E_a(\bar{S}, \bar{\theta}, z)$, which is quite different from Zemskova et al. (2015) or MacCready and Giddings (2016), thus paving the way for improving the description of the Lorenz energy cycle in ocean models, as will be developed and discussed in subsequent work.

By interpreting the mean and eddy fields as the resolved and unresolved components of a numerical ocean model, we use our framework to discuss and understand the energetic implications of different possible parameterisations of the turbulent fluxes affecting the across-scale KE and APE energy transfers, similarly to the energetically consistent modelling approach of Carsten Eden (Eden et al., 2014; Eden, 2015, 2016) (note, though, that in the latter studies, potential energy is partitioned into dynamic and potential enthalpy rather than into APE and BPE). Our main focus is on advective and diffusive parameterisations of the turbulent density flux, which in the oceanic case pertain to meso-scale and small scales, but which in our framework cannot be distinguished. Our discussion reveals that some conclusions about the impact of turbulent parameterisations depend on whether the exact or QG approximation is used; it also suggests that such parameterisations can occasionally convert unresolved APE into resolved APE, potentially

revealing a new type of numerical diffusive instability. We also find that if the energy transfers associated with the APE cascade are assumed to scale proportionally to the KE cascade, aspects of the parameterisations of vertical momentum transfer discussed by Greatbatch and Lamb (1990) can be recovered.

In line with classical discussions of the Lorenz energy cycle, our exact mean/eddy decomposition framework relies on conventional Eulerian Reynolds averaging. However, we acknowledge that there is increasing interest in understanding the energetics for the thickness-weighted averaged equations in isentropic or isopycnic coordinates (e.g., Bleck, 1985; Aiki et al., 2016; Loose et al., 2023). How to extend our results to this case is left to future work.

This paper is organised as follows. Section 2 highlights the convexity of the local APE density as the fundamental property underlying the construction of the concept of eddy APE in the most general case. The convexity property of the local APE density was only briefly mentioned by Scotti and White (2014) but arguably warrants a more thorough discussion and exploitation. This section also clarifies the links between exact and heuristic forms of local APE density. Section 3 revisits and simplifies the derivation of the local budgets of mean and eddy APE previously obtained by Scotti and White (2014). Section 4 discusses the constraints on mixing parameterisations derived from the consideration of the eddy APE and KE budgets, a key issue for the development of energetically consistent numerical ocean models, which do not appear to have been considered before. Section 5 provides a summary and discussion of the results.

## 2. Convexity and eddy APE density

### 2.1. Boussinesq model equations

In this study, we analyse the energetics of rotating stratified flows using the standard Boussinesq approximation. We define the system's state relative to the Lorenz reference state, which represents the configuration of minimum potential energy achievable through an adiabatic rearrangement of fluid parcels. This reference state is characterised by the pressure and density profiles, $p_0(z)$ and $\rho_0(z) = -g^{-1}dp_0/dz$, respectively. With these assumptions, the governing equations of motion may be written as

$$ \frac{D\mathbf{v}}{Dt} + 2\Omega \times \mathbf{v} + \frac{1}{\rho_*} \nabla p_\ell = b_\ell \mathbf{k} + \nu \nabla^2 \mathbf{v} \quad (7) $$

$$ \nabla \cdot \mathbf{v} = 0 \quad (8) $$

$$ \frac{D\rho}{Dt} = -\nabla \cdot \mathbf{J}_{\rho}, \quad \mathbf{J}_{\rho} = -\kappa \nabla \rho. \quad (9) $$
where $p_\ell = p - p_0(z)$ is the pressure anomaly relative to the reference pressure, $v = (u, v, w)$ is the three-dimensional velocity field, $\Omega$ is the rotation vector, $p$ is the pressure, $\rho$ is the density, $\kappa$ is the molecular diffusivity, $v$ is the kinematic viscosity, $\rho_*$ is the constant reference Boussinesq density, and $g$ is the acceleration of gravity.

The buoyancy term in (7) is defined as

$$b_{\ell} = - \frac{g(\rho - \rho_0(z))}{\rho_*} \qquad (10)$$

This buoyancy $b_\ell$ differs from the standard buoyancy $b_{bou} = -g(\rho - \rho_*)/\rho_*$ as it is measured relative to the variable reference density $\rho_0(z)$ rather than a constant Boussinesq density $\rho_*$. For the purpose of this analysis, we assume that the overall domain, analogous to oceanic conditions, is sufficiently large that the reference density profile $\rho_0(z)$ can be considered time-independent. This assumption is supported by climatological observations of temperature and salinity over a century (not shown), which indicate that below the mixed layer, the Lorenz reference state appears to remain stable over time.

## 2.2. Local APE theory

The local APE theory, which builds upon the global APE theory by Lorenz (1955), was initially developed by Andrews (1981) and Holliday and McIntyre (1981), and later rooted in Hamiltonian theory by Shepherd (1993). This theory has been further extended and refined by Scotti et al. (2006), Roullet and Klein (2009), Scotti and White (2014), Zemskova et al. (2015) and Tailleux (2013b, 2018) among others. For a comprehensive review, see Tailleux (2013a). Unlike Lorenz's global APE theory, the local APE theory defines APE as a local, non-negative quantity, expressed through an APE density function whose precise form depends on the equation of state and the approximations used.

In this paper, we specifically consider a standard Boussinesq fluid with a linear equation of state. The expression for the APE density, derived from Holliday and McIntyre (1981) and subsequently utilised by Roullet and Klein (2009) and Tailleux (2013b), is:

$$E_a(\rho, z) = \frac{g}{\rho_*} \int_{z_0(\rho)}^{z} [\rho - \rho_0(\bar{z})] d\bar{z} = - \int_{z_0(\rho)}^{z} b_{\ell}(\rho, \bar{z}) d\bar{z} \qquad (11)$$

where $z_0(\rho)$ is the Level of Neutral Buoyancy (LNB) at which the density of a fluid parcel equals that of the reference density:

$$\rho = \rho_0(z_0(\rho)). \qquad (12)$$

Note that in Eq. (11), $\rho-\rho_0(\bar{z})$ should be interpreted as $\rho(x, y, z, t)-\rho_0(\bar{z})$, with $\rho$ held constant during the integration. The APE density, like other forms of exergy (Marquet, 1991; Kucharski, 1997, 2001), is an extrinsic state function, dependent on both the fluid parcel's state and its environmental context.

The APE density can also be expressed in terms of $\rho$:

$$E_a(\rho, z) = \int_{\rho^2(z)}^{\rho} \frac{\partial E_a}{\partial \rho}(\bar{\rho}, z) d\bar{\rho} = \frac{g}{\rho_*} \int_{\rho^2(z)}^{\rho} [z - z_0(\bar{\rho})] d\bar{\rho}, \qquad (13)$$

where $\rho^\ddagger$ is defined such that $z_0(\rho^\ddagger) = z$. Eq. (13) is the starting point of the mean/eddy decomposition obtained by Scotti and White (2014) and is often favoured over (11) in the literature. In this study, Eq. (11) is preferred over (13), as it aligns more closely with the APE density for a multi-component compressible fluid, e.g., Tailleux (2013b, 2018).

Given that $\rho$ is a function of position and time, $E_a(\rho, z)$ can also be viewed as a function of $(x, y, z, t)$. To differentiate between vertical derivatives calculated at constant $(x, y, t)$ versus at constant $\rho$, we introduce the two separate notations:

$$\frac{\partial}{\partial z} = \frac{\partial}{\partial z}\bigg|_{x,y,t} \quad \text{versus} \quad \frac{\partial}{\partial Z} = \frac{\partial}{\partial z}\bigg|_{\rho}. \qquad (14)$$

The partial derivatives of $E_a$ with respect to density and height are:

$$\frac{\partial E_a}{\partial \rho} = \frac{g(z - z_0(\rho))}{\rho_*} = \frac{g\zeta}{\rho_*} = Y \qquad (15)$$

$$\frac{\partial E_a}{\partial Z} = \frac{g(\rho - \rho_0(z))}{\rho_*} = -b_{\ell}. \qquad (16)$$

Here, $\zeta$ represents the displacement from the reference depth $z_0(\rho)$, and $Y$ denotes a thermodynamic efficiency factor, indicating how diabatic heating influences APE density versus background potential energy (BPE). This concept aligns with the thermodynamic efficiency in compressible fluids, as discussed by Tailleux (2024).

The functional dependence of $E_a$ on $\rho$ and $z$ leads to the following identity

$$\nabla E_a = \frac{\partial E_a}{\partial \rho} \nabla \rho + \frac{\partial E_a}{\partial Z} \nabla z = \mathbf{P}_a - b_{\ell} \mathbf{k} \qquad (17)$$

in which

$$\mathbf{P}_a = \frac{\partial E_a}{\partial \rho} \nabla \rho = Y \nabla \rho. \qquad (18)$$

This equality represents a particular instance of the Crocco-Vazsonyi theorem (Crocco, 1937; Vazsonyi, 1945), which is crucial for developing energetically consistent sound-proof approximations, as recently explored by Tailleux and Dubos (2024). Here, $\mathbf{P}_a$ denotes an APE-based modification of the P-vector originally introduced by Nycander (2011) and recently showed by Tailleux and Wolf (2023) to relate to the directions of lateral stirring in the oceans.

## 2.3. 'Heuristic' eddy APE

In most studies, buoyancy forces are often introduced in an ad-hoc manner, reliant on an arbitrary selection of a reference state (Thorpe et al., 1989; Smith et al., 2005). These forces, while useful, typically lack intrinsic physical significance. However, the work done against the buoyancy forces defined relative to the Lorenz reference density profile represents the energy required to achieve the stratification of the actual state from the Lorenz reference state through an adiabatic rearrangement of fluid parcels. This process has intrinsic dynamical significance since it is inherently tied to the system's physical properties. The associated squared buoyancy frequency profile for such forces is

$$N_0^2(z) = - \frac{g}{\rho_*} \frac{d\rho_0}{dz}(z). \qquad (19)$$

However, it is generally considered that this definition of buoyancy forces pertains primarily to the energetics of the large-scale flows, with less relevance to local turbulent mixing events. For the latter, it is generally assumed that the relevant buoyancy forces are those defined in terms of the buoyancy anomaly:

$$b' = - \frac{g}{\rho_*} (\rho - \bar{\rho}(x, y, z)) = - \frac{g \rho'}{\rho_*}. \qquad (20)$$

Based on this, an intuitive extension of the APE density to measure the work against these small-scale forces consists in adapting (11) as follows:

$$E_a^{\text{heu}} = \frac{g}{\rho_*} \int_{z_m(x, y, \rho)}^{z} [\rho - \bar{\rho}(x, y, \bar{z})] d\bar{z}, \qquad (21)$$

Here, $\bar{\rho}$ replaces $\rho_0(z)$, and $z_m = z_m(x, y, \rho)$ is the height at which a fluid has zero buoyancy relative to $\bar{\rho}$, solving

$$\bar{\rho}[x, y, z_m(x, y, \rho)] = \rho. \qquad (22)$$

This heuristic APE density underpins the work of Roullet et al. (2014) on mesoscale eddy APE from ARGO float data, and Luecke et al. (2017)'s comparison of simulated and observed eddy APE. At leading order, this can be approximated by

$$E_a^{\text{heu}} \approx -\frac{1}{2} b' \zeta' \approx \frac{g^2}{\rho_*^2 \bar{N}^2} \frac{\rho'^2}{2} = \frac{1}{2} \frac{b'^2}{\bar{N}^2} \qquad (23)$$

where $\zeta' = z - z_m(x, y, \rho)$ and $\bar{N}^2$ is the local mean squared buoyancy frequency

$$\bar{N}^2 = - \frac{g}{\rho_*} \frac{\partial \bar{\rho}}{\partial z}(x, y, z_m). \qquad (24)$$
However, defining the APE density of small or mesoscale motions as per (21) has issues:

**Energy budget consistency:** the evolution equation for $E_a^{heu}$ includes energy conversion terms (such as one proportional to the horizontal density gradient $\nabla_h \bar{\rho}$ for instance) with no counterpart in the eddy kinetic energy equation, complicating its energetics description.

**Link to large-scale APE** Physically, the APE dissipation by small scale mixing must ultimately balance, if only partially, the large-scale sources of APE imparted at the boundaries of the domain (Zemskova et al., 2015). Since (21) does not depend on $\rho_0(z)$, it is challenging to link small-scale turbulent dissipation to large scale APE sources.

In this paper, we argue that in order to retain a connection with both $b'$ and $\rho_0(z)$, the eddy APE density needs to be defined as part of an exact mean/eddy decomposition of the APE density. At leading order, this exact eddy APE density is also equal to $-1/2b'\zeta'$, but with the displacement $\zeta'$ defined as $\zeta' = z_0(\rho) - z_0(\bar{\rho})$. Importantly, our exact eddy APE density is valid for arbitrarily large deviations from Lorenz reference state. As established further in the text, this theory modifies the classical relationship between turbulent diapycnal diffusivity and eddy APE dissipation $\epsilon_p^t$:

$$K_{\rho} = \frac{\epsilon_{\rho}^{t}}{N^2} \quad (old) \quad \Rightarrow \quad K_{\rho} = \frac{1}{\Lambda(1 + |S|^2)} \frac{\epsilon_{\rho}^{t}}{N^2} \quad (new) \quad (25)$$

Here, $\Lambda$ and $S$ are parameters defined by

$$\Lambda = \frac{\partial \bar{\rho}}{\partial z} \frac{\partial z_0}{\partial \rho} (\bar{\rho}) = \frac{\bar{N}^2}{\bar{N}_0^2}, \quad S = - \left( \frac{\partial \bar{\rho}}{\partial z} \right)^{-1} \nabla_h \bar{\rho} \quad (26)$$

where $\bar{N}_0^2 = N_0^2(z_0(\bar{\rho}))$. The two parameters $\Lambda$ and $S$ measure deviations from Lorenz reference state in different ways and play a crucial role in this paper, as clarified further in the text. Neglecting $|S|$ in the expression for $K_\rho$ is equivalent to making the small slope approximation. Physically, $\Lambda$ is expected to differ significantly from unity where the local squared buoyancy differs significantly from its value in Lorenz reference state; in the oceans, this is primarily the case in the polar oceans (Saenz et al., 2015; Tailleux, 2016; Tailleux and Wolf, 2023).

## 2.4. Convexity of APE density and eddy APE

We now separate variables into mean and eddy components using standard Reynolds averaging so that $v = \bar{v} + v'$, $\rho = \bar{\rho} + \rho'$, and so on, such that for any variable $f$, $\bar{f} = \bar{f}'$ and $f' = 0$. As is well known, such an approach yields to the following mean/eddy partition of the Reynolds averaged kinetic energy

$$\frac{\bar{v}^2}{2} = \frac{\bar{v}^2}{2} + \frac{\bar{v}^2}{2} = E_k^m + E_k^t. \quad (27)$$

Importantly, the mean kinetic energy $E_k^m$ appears as the kinetic energy of the mean velocity field $\bar{v}$. The corresponding problem for the APE density is to achieve the following mean/eddy decomposition of the Reynolds averaged APE density

$$\bar{E}_a = E_a^m + E_a^t. \quad (28)$$

Ideally, one would like to define the mean APE density as $E_a^m = E_a(\bar{\rho}, z)$, that is, as the APE density of the mean density $\bar{\rho}$, since the latter is the quantity that appears in the Reynolds averaged momentum equations, and therefore define the eddy APE as the residual $E_a^t = \bar{E}_a - E_a(\bar{\rho}, z)$. That such an approach leads to a non-negative $E_a^t$ was established by Scotti and White (2014), who linked the result to the convexity of APE density with respect to buoyancy but without elaborating on it. However, given its fundamental importance in available energy theories, the convexity property warrants more emphasis.

Mathematically, a function $f(x)$ is said to be convex at some point $x_0$ if its curve lies above its tangent line at that point, hence if the quantity $f_e(x; x_0) = f(x) - f(x_0) - f'(x_0)(x - x_0) \ge 0$ is non-negative. To clarify what properties of $f$ determine its convexity, it is useful to rewrite $f_e$ in the form

$$f_e = \int_{x_0}^{x} [f'(\bar{x}) - f'(x_0)] d\bar{x} = \int_{x_0}^{x} \int_{x_0}^{\bar{x}} f''(\bar{x}) d\bar{x} d\bar{x}. \quad (29)$$

Eq. (29) is an important identity as it shows that $f$ is convex at the point $x_0$ if its second derivative $f''$ is non-negative (assuming $f$ to be twice differentiable). Of course, convexity extends to functions of several variables. For instance, for a function $f(x, y)$, convexity requires that the quantity $f_e(x, y; x_0, y_0) = f(x, y) - f(x_0, y_0) - \partial_x f(x_0, y_0)(x - x_0) - \partial_y f(x_0, y_0)(y - y_0) \ge 0$ be non-negative. Convexity plays a key role in thermodynamics. For instance, it can be shown that the possibility to convert heat into work, the central object of thermodynamics, hinges on internal energy being a convex function of its canonical variables, specific entropy and specific volume, which is key to defining the concept of exergy, e.g., Tailleux (2013a). Arguably, it is the convexity of kinetic energy that is implicitly responsible for the non-negative character of eddy kinetic energy. For quadratic expressions, however, there is no need to invoke convexity as this is not required to prove the non-negative character of the eddy component. Convexity is needed here, however, because the APE density includes higher-order non-quadratic terms, called anharmonic by Roullet and Klein (2009).

As it turns out, the APE density (11) is convex with respect to both density and $z$, which can be verified by differentiating (15) and (16) with respect to $\rho$ and $z$ respectively, which leads to

$$\frac{\partial^2 E_a}{\partial \rho^2} = - \frac{g}{\rho_*} \frac{d z_0}{d \rho}(\rho), \quad (30)$$

$$\frac{\partial^2 E_a}{\partial Z^2} = - \frac{g}{\rho_*} \frac{d \rho_0}{d z}(z) = N_0^2(z). \quad (31)$$

Physically, the non-negative character of $N_0^2$ in Eq. (16) follows from Lorenz reference state being a state of minimum potential energy, hence statically stable by construction, which establishes convexity with respect to $z$. To prove the non-negative character of (15) and hence the convexity with respect to density, simply differentiate the LNB Eq. (12) with respect to $\rho$, thus leading to

$$\frac{d \rho_0}{d z}(z_0(\rho)) \frac{d z_0}{d \rho}(\rho) = 1. \quad (32)$$

Now, since $d\rho_0/dz \le 0$, (32) implies that $dz_0/d\rho \le 0$, which proves our proposition.

Having established the convexity of $E_a$ with respect to both density and $z$, let us define the instantaneous eddy APE as the non-negative nonlinear term $A_e$ in perturbation density in the following series expansion of $E_a$ around $\bar{\rho}$

$$E_a(\rho, z) = E_a(\bar{\rho}, z) + \frac{\partial E_a}{\partial \rho}(\bar{\rho}, z)(\rho - \bar{\rho}) + A_e. \quad (33)$$

Upon Reynolds averaging, the term proportional to $\rho' = \rho - \bar{\rho}$ vanishes, leading the mean/eddy decomposition:

$$\bar{E}_a = E_a^m + E_a^t, \quad (34)$$

with

$$E_a^m = E_a(\bar{\rho}, z), \quad E_a^t = \bar{A}_e. \quad (35)$$

Physically, the instantaneous eddy APE $A_e$ is the counterpart of the instantaneous value of eddy kinetic energy $v^2/2$. To facilitate the
comparison of $A_e$ with the heuristic APE density discussed earlier, it is useful to rewrite $A_e$ in the following more revealing form

$$ 
\begin{aligned}
A_e &= E_a(\rho, z) - E_a(\bar{\rho}, z) - \frac{\partial E_a}{\partial \bar{\rho}}(\bar{\rho}, z)(\rho - \bar{\rho}) \\
&= \int_{\bar{\rho}}^{\rho} \left[ \frac{\partial E_a}{\partial \rho}(\bar{\rho}, z) - \frac{\partial E_a}{\partial \rho}(\bar{\rho}, z) \right] d\bar{\rho} = - \frac{g}{\rho_*} \int_{\bar{\rho}}^{\rho} [z_0(\bar{\rho}) - z_0(\rho)] d\bar{\rho}
\end{aligned}
\quad (36) $$

where the passage from the penultimate to last equation made use of (15). Alternatively, $A_e$ may be rewritten as an integral of the work against buoyancy forces by introducing the change of variable $\bar{\rho} = \rho_0(\bar{z})$ so that $d\bar{\rho} = \rho'_0(\bar{z})d\bar{z}$, and

$$ A_e = - \frac{g}{\rho_*} \int_{z_0(\bar{\rho})}^{z_0(\rho)} [\rho - \rho_0(\bar{z})] d\bar{z} \quad (37) $$

Importantly, note that (36) and (37) reveal that $A_e = A_e(\rho, \bar{\rho})$ no longer depends on height $z$, being solely a function of $\rho$ and $\bar{\rho}$ only. In both cases, the expression for the transient eddy APE may be approximated as

$$ A_e \approx -\frac{1}{2} b' \zeta' \quad \Rightarrow \quad E_a' \approx -\frac{1}{2} \bar{b}' \zeta' \quad (38) $$

with

$$ b' = - \frac{g(\rho - \bar{\rho})}{\rho_*}, \quad \zeta' = z_0(\rho) - z_0(\bar{\rho}) \quad (39) $$

The approximation (38) is identical to the expression used by Roullet et al. (2014). However, while the buoyancy anomaly $b'$ is the same as in Roulet's approach, the displacement $\zeta'$ is defined quite differently in terms of the instantaneous and mean reference positions $z_0(\rho)$ and $z_0(\bar{\rho})$ respectively, that is, in terms of the equilibrium positions of $\rho$ and $\bar{\rho}$ in Lorenz reference state. These differences are clearly illustrated in the right panel of Fig. 1 introduced earlier.

Physically, the approximation (38) is obtained from a simple trapezoidal approximation of the integrals (36) or (37), and are likely to be the most accurate approximation of $A_e$. Nevertheless, using the approximation

$$ \zeta' = z_0(\rho) - z_0(\bar{\rho}) \approx \frac{\partial z_0}{\partial \rho}(\bar{\rho})\rho' \quad (40) $$

it is also possible to approximate $A_e$ in terms of the following quadratic expressions:

$$ A_e \approx -\frac{g}{\rho_*} \frac{dz_0}{d\rho}(\bar{\rho}) \frac{\rho'^2}{2} = \frac{g^2}{\rho_*^2 \bar{N}_0^2} \frac{\rho'^2}{2} = \frac{1}{2} \frac{b'^2}{\bar{N}_0^2} \quad (41) $$

with

$$ \bar{N}_0^2 = \bar{N}_0^2(x, y, z) = - \frac{g}{\rho_*} \frac{d\rho_0}{dz}(z_0(\bar{\rho})). \quad (42) $$

Note that $\bar{N}_0^2$, unlike $\bar{N}_0^2(z)$, is a function of all three spatial dimensions $(x, y, z)$. Its spatial gradient is easily verified to be

$$ \nabla \bar{N}_0^2 = - \frac{g}{\rho_*} \frac{d^2 \rho_0}{dz^2} \frac{\partial z_0}{\partial \rho} \nabla \bar{\rho} \quad (43) $$

so appears to be proportional to the mean density gradient $\nabla \bar{\rho}$, the proportionality factor being controlled by the curvature of $\rho_0(z)$. Comparison with the heuristic localised APE density $E_a^{heu}$ and $A_e$ is easily verified to be

$$ \frac{A_e}{E_a^{heu}} \approx \frac{\bar{N}_0^2}{\bar{N}_0^2} = \Lambda, \quad (44) $$

where $\Lambda$ is the same parameter introduced previously. This establishes that the validity and accuracy of the heuristic localised APE density depend on the proximity of the actual state to Lorenz reference state, as expected.

### 3. Local budgets of available potential energy

#### 3.1. Non-averaged local APE budgets

Before deriving local budget equations for the mean and APE densities $E_a^m$ and $E_a'$, we first clarify the local budget equation satisfied by the non-averaged APE density $E_a$. This can be obtained by taking the Lagrangian derivative of $E_a$, yielding

$$ \frac{DE_a}{Dt} = \left( \frac{\partial E_a}{\partial \rho} \frac{D\rho}{Dt} + \frac{\partial E_a}{\partial Z} \frac{Dz}{Dt} \right) = Y \frac{D\rho}{Dt} - b_e w. \quad (45) $$

By making use of the density Eq. (9), it is easily checked that (45) may be rewritten in the form

$$ \frac{DE_a}{Dt} = -b_e w - \nabla \cdot \mathbf{J}_a - \varepsilon_p \quad (46) $$

in which $\mathbf{J}_a$ and $\varepsilon_p$ are the diffusive flux of APE density and APE dissipation rate, respectively, given by

$$ \mathbf{J}_a = Y \mathbf{J}_p = -Y \kappa \nabla \rho, \quad (47) $$

$$ \varepsilon_p = -\mathbf{J}_p \cdot \nabla Y = \kappa \nabla \rho \cdot \nabla Y \quad (48) $$

The normal component of $\mathbf{J}_a$ at the ocean surface determines the APE production rate by surface buoyancy fluxes, see Zemskova et al. (2015) for a discussion within the present framework. Our evolution equation for the APE density (46), although mathematically equivalent, is much simpler in form that the one previously derived (Scotti and White, 2014) due to not imposing $\mathbf{J}_a$ to be downgradient in $E_a$. Physically, the form (46) is to be preferred because it is the one that most naturally generalise to double diffusive multi-component compressible stratified fluids, unlike Scotti and White (2014)'s approach, e.g. see Tailleux (2024) for details. Using the expression (15) for $Y$, $\varepsilon_p$ may also be expressed in the more familiar form

$$ \varepsilon_p = \frac{g\kappa}{\rho_*} \nabla \rho \cdot \left[ \mathbf{k} - \frac{dz_0}{d\rho} \nabla \rho \right] = \frac{g\kappa}{\rho_*} \left[ \frac{\partial \rho}{\partial z} - \frac{dz_0}{d\rho} |\nabla \rho|^2 \right] \quad (49) $$

How the APE density budget Eq. (46) relates to that of a fully compressible fluid as well as to the global evolution equations derived by Winters et al. (1995), is extensively discussed by Tailleux (2024), to which the reader is referred to for details. For the reader more familiar with the global APE approach of Winters et al. (1995), it may be useful to point out that the volume integral of (49) coincides with the term $\Phi_d - \Phi_i$ of Winters et al. (1995).

#### 3.2. Mean APE budget

As established previously, the mean APE density $E_a^m = E_a(\bar{\rho}, z)$ is naturally defined as the APE density of the mean density $\bar{\rho}$ and $z$, at least when approached from the viewpoint of standard Eulerian averaging. An evolution equation for it can therefore be derived essentially as that for the non-averaged density $E_a(\rho, t)$, with $D\rho/Dt$ replaced by $D_m\bar{\rho}/Dt$, the Lagrangian derivative of $\bar{\rho}$ defined in terms of the mean velocity $\bar{v}$, which leads to

$$ \frac{D_m E_a^m}{Dt} = \frac{\partial E_a^m}{\partial \bar{\rho}} \frac{D_m \bar{\rho}}{Dt} + \frac{\partial E_a^m}{\partial Z} \frac{D_m z}{Dt} = Y_m \frac{D_m \bar{\rho}}{Dt} - \bar{b}_e \bar{w}, \quad (50) $$
in which $Y_m$ and $\bar{b}_\ell$ are defined by

$$Y_m = \frac{\partial E_a}{\partial \rho} (\bar{\rho}, z) = \frac{g(z - z_0(\bar{\rho}))}{\rho_*} \qquad (51)$$

$$\bar{b}_\ell = \frac{\partial E_a}{\partial Z}(\bar{\rho}, z) = \frac{g(\bar{\rho} - \rho_0(z))}{\rho_*} \qquad (52)$$

It is important to note here that the mean quantity $Y_m$ differs from the Reynolds averaged $Y$, i.e., $Y_m \neq \bar{Y}$, because $z_0(\bar{\rho}) \neq z_0(\rho)$ in general. Using the evolution equation for the mean density $\bar{\rho}$, viz.,

$$\frac{D_m \bar{\rho}}{Dt} = -\nabla \cdot [\overline{\rho' \mathbf{v}'} - \kappa \nabla \bar{\rho}] = -\nabla \cdot \mathbf{J}_\rho^m \qquad (53)$$

in which $\mathbf{J}_\rho^m = \overline{\rho' \mathbf{v}'} - \kappa \nabla \bar{\rho}$ is the total density flux including both turbulent and molecular diffusive contributions, it is straightforward to show that (50) may be rewritten as

$$\frac{D_m E_a^m}{Dt} = -\bar{b}_\ell \overline{\mathbf{w}} - \nabla \cdot \mathbf{J}_a^m + \mathbf{J}_\rho^m \cdot \nabla Y_m \qquad (54)$$

where $\mathbf{J}_a^m$ is the total ‘diffusive’ flux of mean APE density including both turbulent and molecular contributions, given by

$$\mathbf{J}_a^m = Y_m \mathbf{J}_\rho^m = Y_m (\overline{\rho' \mathbf{v}'} - \kappa \nabla \bar{\rho}) \qquad (55)$$

Note that the last term in (54) may be further expanded in the form

$$\mathbf{J}_\rho^m \cdot \nabla Y_m = \overline{\rho' \mathbf{v}'} \cdot \nabla Y_m - \kappa \nabla \bar{\rho} \cdot \nabla Y_m = \overline{\rho' \mathbf{v}'} \cdot \nabla Y_m - \epsilon_\rho^m \qquad (56)$$

where

$$\epsilon_\rho^m = \kappa \nabla \bar{\rho} \cdot \nabla Y_m \qquad (57)$$

represents the ‘mean’ APE dissipation rate. As a result, (54) may ultimately be rewritten in the following final form:

$$\frac{D_m E_a^m}{Dt} = -\bar{b}_\ell \overline{\mathbf{w}} + \overline{\rho' \mathbf{v}'} \cdot \nabla Y_m - \nabla \cdot \mathbf{J}_a^m - \epsilon_\rho^m \qquad (58)$$

Physically, the terms appearing in the r.h.s. of (58) represent: (1) the conversion between mean APE and mean KE; (2) the conversion between mean APE and eddy APE; (3) the diffusive flux of mean APE by means of turbulent and molecular processes; (4) the mean dissipation rate of APE by molecular processes.

### 3.3. Eddy APE budget and turbulent APE dissipation

We now turn to the problem of deriving an evolution equation for the eddy APE $E_a^t = \overline{E}_a$. There are two main routes. The first route is via deriving an evolution equation for $A_e$ and Reynolds averaging the result. In the second route, which is much simpler and the only one pursued here, the evolution equation for $E_a^t = \overline{E}_a - E_a^m$ is obtained as the residual between the evolution equations for $\overline{E}_a$ and $E_a^m$.

To proceed, we first take the Reynolds average of (46) after separating each variable into mean and eddy components, which leads to

$$\frac{D_m \overline{E}_a}{Dt} + \nabla \cdot (\overline{E_a' \mathbf{v}'}') = -\overline{b}_\ell \overline{\mathbf{w}} - \overline{b'_\ell \mathbf{w}'} - \nabla \cdot \overline{\mathbf{J}_a} - \overline{\epsilon}_\rho \qquad (59)$$

The sought-for evolution equation for $E_a^t$ is then simply obtained by subtracting the $E_a^m$ Eq. (58) from (59), which yields

$$\frac{D_m E_a^t}{Dt} = -\overline{b'_\ell \mathbf{w}'} - \overline{\rho' \mathbf{v}'} \cdot \nabla Y_m - \nabla \cdot \mathbf{J}_a^t - \epsilon_\rho^t \qquad (60)$$

in which $\mathbf{J}_a^t$ and $\epsilon_\rho^t$ are the total flux of eddy APE density and eddy APE dissipation rate respectively, whose expressions are

$$\mathbf{J}_a^t = \overline{\mathbf{J}_a} - \mathbf{J}_a^m + \overline{E_a' \mathbf{v}'} \qquad (61)$$

$$\epsilon_\rho^t = \overline{\epsilon}_\rho - \epsilon_\rho^m \qquad (62)$$

Using the fact that $\mathbf{J}_a^m = -Y_m \kappa \nabla \bar{\rho}$, $\epsilon_\rho^m = \kappa \nabla \bar{\rho} \cdot \nabla Y_m$, and $\overline{\epsilon}_\rho = \kappa \nabla \bar{\rho} \cdot \overline{\nabla Y} + \kappa \nabla \bar{\rho}' \cdot \overline{\nabla Y}'$, these can be more explicitly written as

$$\mathbf{J}_a^t = -\kappa \overline{Y' \nabla \rho'} - (\overline{Y} - Y_m) \kappa \nabla \bar{\rho} + \overline{E_a' \mathbf{v}'} \qquad (63)$$

$$\epsilon_\rho^t = \overline{\epsilon}_\rho - \epsilon_\rho^m = \kappa \nabla \bar{\rho} \cdot \nabla (\overline{Y} - Y^m) + \kappa \overline{\nabla \rho'} \cdot \nabla Y' \qquad (64)$$

These expressions show that the turbulent diffusive flux of APE and APE dissipation rate both depend on the mean quantity $\overline{Y} - Y_m$, whose leading order expression can be showed to be

$$\overline{Y} - Y_m = \frac{g}{\rho_*} (\overline{z_0(\rho)} - z_0(\bar{\rho})) \approx \frac{g}{\rho_*} \frac{\partial^2 z_0}{\partial \rho^2} (\bar{\rho}) \frac{\rho'^2}{2} \qquad (65)$$

The above expressions also depends on $Y'$, which can be shown to be given at leading order

$$Y' = \frac{g}{\rho_*} (\overline{z_0(\rho)} - z_0(\bar{\rho})) \approx -\frac{g}{\rho_*} \frac{\partial z_0}{\partial \rho} (\bar{\rho}) \rho' + \frac{g}{\rho_*} \frac{\partial^2 z_0}{\partial \rho^2} (\bar{\rho}) \left(\frac{\overline{\rho'^2}}{2} - \frac{\rho'^2}{2}\right) \qquad (66)$$

As a result, it follows that

$$\kappa \nabla \bar{\rho} \cdot \nabla (\overline{Y} - Y_m) \approx \frac{g}{\rho_*} \frac{\partial^2 z_0}{\partial \rho^2} (\bar{\rho}) \kappa \nabla \bar{\rho} \cdot \nabla \frac{\overline{\rho'^2}}{2} + \frac{g}{\rho_*} \frac{\overline{\rho'^2}}{2} \frac{\partial^3 z_0}{\partial \rho^3} (\bar{\rho}) \kappa |\nabla \bar{\rho}|^2 \qquad (67)$$

$$\kappa \overline{\nabla \rho'} \cdot \nabla Y' \approx - \frac{g}{\rho_*} \frac{\partial z_0}{\partial \rho} (\bar{\rho}) \kappa |\overline{\nabla \rho'}|^2 - \frac{g}{\rho_*} \frac{\partial^2 z_0}{\partial \rho^2} (\bar{\rho}) \kappa \nabla \bar{\rho} \cdot \nabla \frac{\overline{\rho'^2}}{2} \qquad (68)$$

Summing up these two results, retaining only the terms up to second order in density perturbation, gets rid of the term proportional to $\partial^2 z_0 / \partial \rho^2$ and yields the following equation for the turbulent APE dissipation rate

$$\epsilon_\rho^t \approx -\frac{g}{\rho_*} \frac{\partial z_0}{\partial \rho} (\bar{\rho}) \kappa |\overline{\nabla \rho'}|^2 + \frac{g}{\rho_*} \frac{\overline{\rho'^2}}{2} \frac{\partial^3 z_0}{\partial \rho^3} (\bar{\rho}) \kappa |\overline{\nabla \bar{\rho}}|^2 \qquad (69)$$

Of the two terms appearing in the right-hand side of (69), only the first one can be ascertained to be non-negative and directly comparable to the expression for the APE dissipation rate proportional to the turbulent dissipation rate of density variance previously derived by Oakey (1982) and Gargett and Holloway (1984). The second term, however, can be of any sign, but is expected to be much smaller than the first term in general and therefore unlikely to be important in practice, although this remains to be checked more systematically in direct numerical simulations of turbulent stratified mixing. This is left for future work. Scotti and White (2014) only retained the first term in the right-hand side of (69) in their paper.

## 4. Application to energetically consistent numerical ocean modelling

Our framework provides a rigorous theoretical foundation for linking parameterised energy transfers to observable KE and APE dissipation rates. In numerical ocean models, the main energy transfers of interest are those associated with the ocean meso-scale, which are responsible for transferring mean APE to turbulent or eddy kinetic energy (which is subsequently dissipated through irreversible viscous processes), and those associated with small-scale turbulent stratified mixing, which transfer turbulent kinetic energy into turbulent APE, later dissipated by irreversible diffusive processes, giving rise to turbulent diapycnal mixing. The following aims to illustrate these ideas.

### 4.1. Energetically consistent modelling

Parameterisations of subgridscale processes control the energy transfers between the resolved and unresolved scales of motion implicated in the turbulent forward KE and APE energy cascades. Upon reaching molecular scales, mechanical energy can be dissipated quasi-adiabatically as KE at the viscous dissipation rate $\epsilon_k^t$, or diabatically
as APE at the diffusive dissipation rate $\epsilon_p^t$. In most numerical ocean models, subgrid-scale process parameterisations have traditionally been developed without explicit consideration of their implications for the KE and APE cascades or for the KE and APE dissipation rates. However, this changed when Munk and Wunsch (1998) noted that the intensity of mixing depends on the intensity of the mechanical sources of stirring. In particular, they raised the issue of whether tuning the value of the turbulent diapycnal or vertical mixing coefficients $K_\rho$ or $K_\nu$ to minimise mismatch with observed temperature and salinity fields was necessarily consistent with our understanding of how turbulent mixing is mechanically sustained. To that end, they considered the gravitational potential energy (GPE) budget, assuming a balance between the rate of GPE increase due to mixing at the rate $K_\nu N^2$ and the rate of GPE loss due to high-latitude cooling (neglecting the role of nonlinearities in the equation of state). To link $K_\nu$ to the mechanical sources of stirring, they used the relation $K_\nu N^2 = \Gamma \epsilon_K$ with $\Gamma = 0.2$, and balanced the total viscous dissipation rate with the power input due to winds and tides, dismissing the role of surface buoyancy fluxes (see Tailleux (2010) for a discussion). The work of Munk and Wunsch (1998) proved very influential, prompting the development of ‘energetically consistent modelling’ pioneered by Carsten Eden and his group (Eden et al., 2014; Eden, 2015, 2016), whose aim is to link the intensity of mixing to the sources of stirring by exploiting empirical knowledge about the dissipation ratio $\Gamma = \epsilon_p/\epsilon_k$. In our view, this aim can be more generally defined as seeking to exploit existing and developing knowledge about the inter-relations between the KE and APE budgets, for instance by developing and using a turbulent APE budget in addition to the widely used turbulent KE budget.

In practice, the prevailing approach has primarily relied on using a turbulent kinetic energy (TKE) equation to predict the turbulent viscous dissipation $\epsilon_k^t$, from which the APE dissipation rate $\epsilon_p^t$ and turbulent vertical mixing coefficient $K_\rho = \epsilon_p^t/N^2$ can be inferred, provided that the value of the dissipation ratio $\Gamma$ can also be predicted in some way. For a review of our current understanding of $\Gamma$, see Gregg (2021). An alternative, yet to be developed in oceanography, would be to predict $\epsilon_p^t$ directly from a turbulent APE equation, as proposed in the context of atmospheric boundary layer research by Zilitinkevich and collaborators (e.g., Zilitinkevich et al. (2013)). An important advantage of the eddy APE budget is that it constrains the full turbulent density flux $\overline{\rho'v'}$, whereas the TKE budget depends only on the vertical component $\overline{\rho'w'}$, which is why we think it should receive more attention.

To analyse across-scale energy transfers, a filtering approach is needed to isolate the different scales of interest. In the context of large-scale ocean modelling, Eden (2015) proposed that the subgrid-scale energy should be divided into subreservoirs for meso-scale eddies, internal gravity waves, and turbulent incoherent motions. Here, we only consider a subdivision of energy into mean (resolved) and eddy (unresolved) scales due to the inherent limitations of standard Reynolds averaging. With this in mind, we return to the eddy APE budget (46), which we rewrite as

$$ \frac{\partial E_a^t}{\partial t} + \nabla \cdot (\overline{\nabla} E_a^t + \mathbf{J}_a^t) = C(E_k^t, E_a^t) + C(E_a^m, E_a^t) - \epsilon_p^t \quad (70) $$

where $C(E_k^t, E_a^t)$ and $C(E_a^m, E_a^t)$ are the conversions of turbulent kinetic energy and mean APE into eddy APE, respectively, with

$$ C(E_k^t, E_a^t) = -b'_e w' = \frac{g}{\rho_*} \overline{\rho' w'} \quad (71) $$

$$ C(E_a^m, E_a^t) = -\overline{\rho'v'} \cdot \nabla Y_m = -\frac{g}{\rho_*} \overline{\rho'w'} + \frac{g}{\rho_*} \frac{\partial z_0}{\partial \rho} \overline{(\rho)\overline{\rho'v'}} \cdot \nabla \overline{\rho} \quad (72) $$

Under the classical assumptions of stationarity and homogeneity, the terms on the left-hand side of (70) can be neglected, and the eddy APE budget reduces to a balance between production of eddy APE due to conversions with mean APE and eddy KE, and turbulent APE dissipation, viz.,

$$ C(E_k^t, E_a^t) + C(E_a^m, E_a^t) \approx \epsilon_p^t \quad (73) $$

which, from (71) and (72), may be written as

$$ \frac{g}{\rho_*} \frac{\partial z_0}{\partial \rho} (\overline{\rho} \overline{\rho'v'}) \cdot \nabla \overline{\rho} \approx \epsilon_p^t \quad (74) $$

Eq. (74) can be understood as a constraint on the diapycnal component of the turbulent density flux across the mean isopycnal surfaces $\overline{\rho} = \text{constant}$, controlled by the turbulent APE dissipation rate $\epsilon_p^t$. When the resolved flow pertains to large scales, the turbulent density flux is generally assumed to contain at least two components:

$$ \overline{\rho'v'} = \overline{\rho'v'}_{\text{meso}} + \overline{\rho'v'}_{\text{small}} \quad (75) $$

pertaining to the effects of meso-scale eddies and small-scale turbulent mixing, respectively, and are taken to be perpendicular and parallel to $\nabla \overline{\rho}$ according to

$$ \overline{\rho'v'}_{\text{meso}} = \mathbf{\Psi} \times \nabla \overline{\rho}, \quad \overline{\rho'v'}_{\text{mix}} = -K_\rho \nabla \overline{\rho} \quad (76) $$

(e.g., Griffies, 1998; Griffies et al., 1998). Of particular interest is the vertical component

$$ \overline{\rho'w'}_{\text{mix}} = -K_\rho \frac{\partial \overline{\rho}}{\partial z} \quad (77) $$

which will be later contrasted with the vertical component of the skew diffusive flux.

## 4.2. Energetics of downgradient diffusion

In the literature, the turbulent diffusivity $K_\rho$ entering the turbulent mixing parameterisation for the diffusive part of the turbulent density flux is traditionally predicted by

$$ K_\rho \approx \frac{\epsilon_p^t}{N^2} = \frac{\Gamma \epsilon_k^t}{N^2} \quad (78) $$

(e.g., Lindborg and Brethouwer, 2008), where $\Gamma$ is the dissipation ratio (Oakey, 1982). However, according to the eddy APE budget (76), a more accurate expression is

$$ -\frac{g}{\rho_*} \frac{\partial z_0}{\partial \rho} (\overline{\rho}) K_\rho |\nabla \overline{\rho}|^2 \approx \epsilon_p^t \quad (79) $$

which may be rearranged as

$$ K_\rho \approx \frac{1}{\Lambda(1 + |\mathbf{S}|^2)} \frac{\epsilon_p^t}{\overline{N}^2} \neq \frac{\epsilon_p^t}{\overline{N}^2} \quad (80) $$

where $\Lambda$ and $\mathbf{S}$ are defined as before by (26). As stated earlier, Eq. (80) shows that the standard expression (78) implicitly depends on two assumptions that are rarely, if ever, acknowledged: (1) that the local stratification as measured by $\overline{N}^2$ is approximately equal to $N_0^2$, and (2) the small slope approximation $|\mathbf{S}| \ll 1$, often made in the context of rotated Redi diffusion (Redi, 1982), for instance. Moreover, if we neglect the higher order terms in the definition of the eddy APE dissipation rate so that

$$ \epsilon_p^t \approx -\frac{g}{\rho_*} \frac{\partial z_0}{\partial \rho} (\overline{\rho}) \kappa |\nabla \overline{\rho'}|^2 \quad (81) $$

it is easily verified that the eddy APE budget also implies

$$ K_\rho |\nabla \overline{\rho}|^2 \approx \kappa |\nabla \overline{\rho'}|^2. \quad (82) $$

Physically, Eq. (82) states that the resolved dissipation of the mean density field must ultimately be balanced by the dissipation of eddy density variance at molecular scales. Note, however, that when approached from the eddy APE budget viewpoint, (82) requires neglecting the higher order terms in the eddy APE dissipation rate. In the literature, (82) is more commonly obtained by equating the production of
density variance by the turbulent density flux with the dissipation of density variance, as in the Osborn-Cox model (Osborn and Cox, 1972). Note that Eq. (82) cannot give rise to upgradient (negative) diffusion, while (80) can potentially allow it in some circumstances.

To synthesise the results, the energy conversions affected by the downgradient part of the turbulent density flux are therefore given by

$$ \{C(E_a^m, E_a^t)\}_{\text{mix}} = ((1 + |S|^2)\Lambda - 1) K_\rho \overline{N}^2 \quad (83) $$

$$ \{C(E_k^t, E_a^t)\}_{\text{mix}} = K_\rho \overline{N}^2. \quad (84) $$

In the literature, the mean to eddy APE conversion has been exclusively discussed in the context of the quasi-geostrophic (QG) approximation (e.g., von Storch et al., 2012), which is equivalent to assuming $\Lambda \approx 1$ in (83), in which case it reduces to

$$ \{C(E_a^m, E_a^t)\}_{\text{qg,mix}} \approx |S|^2 K_\rho \overline{N}^2 \ge 0, \quad (85) $$

and always acts as a downscale transfer of energy from mean APE to eddy APE. The exact finite-amplitude mean APE to eddy APE conversion may occasionally behave quite differently, as Eq. (83) shows that the sign of the conversion is no longer necessarily non-negative. Energy transfer can be either upscale or downscale, depending on the relative values of $\Lambda$ and the slope parameter $|S|$. Because the value of $\Lambda$ is controlled by the distance from the Lorenz reference state, which itself depends on the steepness of isopycnal slopes, we assume that the condition for upscale transfer from eddy APE to mean APE imposes a constraint on $\Lambda$ (rather than on $|S|$), namely

$$ \Lambda < \frac{1}{1 + |S|^2}. \quad (86) $$

In that case, downgradient diffusion can, at least in principle, backscatter unresolved energy into resolved energy, thus potentially acting as a source of instability for the resolved flow. Although the associated upscale energy transfer will be counteracted by the downscale energy transfer associated with skew diffusion, it is not necessarily obvious that this can actually suppress the diffusive instability, because downgradient diffusion and skew diffusion a priori operate differently.

### 4.3. Energetics of skew diffusion and TKE budget

Physically, skew diffusion can be interpreted as an eddy-induced advection by subgrid-scale processes, as follows from the relation

$$ \nabla \cdot (\boldsymbol{\Psi} \times \nabla \bar{\rho}) = (\nabla \times \boldsymbol{\Psi}) \cdot \nabla \bar{\rho} = \mathbf{v}_{\text{eddy}} \cdot \nabla \bar{\rho}. \quad (87) $$

As a result, the evolution equation for the mean density $\bar{\rho}$ may be written as

$$ \frac{\partial \bar{\rho}}{\partial t} + \mathbf{v}_{\text{res}} \cdot \nabla \bar{\rho} = \nabla \cdot [(K_{\rho} + \kappa) \nabla \bar{\rho}] \quad (88) $$

with $\mathbf{v}_{\text{res}} = \bar{\mathbf{v}} + \mathbf{v}_{\text{eddy}}$ being the residual velocity and $\mathbf{v}_{\text{eddy}} = \nabla \times \boldsymbol{\Psi}$ the eddy-induced velocity. By construction, the eddy-induced velocity is divergenceless, $\nabla \cdot \mathbf{v}_{\text{eddy}} = 0$. In oceanography, the most commonly used parameterisation for $\boldsymbol{\Psi}$ is

$$ \boldsymbol{\Psi} = \boldsymbol{\Psi}_{\text{gm}} = \mathbf{k} \times \kappa_{\text{gm}} \mathbf{S}, \quad (89) $$

as originally proposed by Gent and McWilliams (1990), Gent et al. (1995), where $\mathbf{S}$ is the slope vector previously introduced. For the GM parameterisation, the skew-diffusion part of the turbulent density flux becomes

$$ \overline{\rho' \mathbf{v}'}_{\text{skew}} = (\mathbf{k} \times \kappa_{\text{gm}} \mathbf{S}) \times \nabla \bar{\rho} = \kappa_{\text{gm}} \frac{\partial \bar{\rho}}{\partial z} \mathbf{S} + \kappa_{\text{gm}} |S|^2 \frac{\partial \bar{\rho}}{\partial z} \mathbf{k}. \quad (90) $$

Of particular interest is the vertical component, given by

$$ \overline{\rho' \omega'}_{\text{skew}} = \kappa_{\text{gm}} |S|^2 \frac{\partial \bar{\rho}}{\partial z}. \quad (91) $$

which acts as an ‘upgradient’ flux.

The skew-diffusive part of the turbulent density flux is perpendicular to the mean density gradient $\nabla \bar{\rho}$, and therefore does not have a net contribution to the eddy APE budget. However, it is associated with net energy conversions between mean and eddy APE, as well as between eddy APE and eddy KE. From (71) and (91), it is easily verified that

$$ \{C(E_a^m, E_a^t)\}_{\text{skew}} = -\overline{\rho' \mathbf{v}'}_{\text{skew}} \cdot \nabla Y_m = -\frac{g}{\rho_*} \overline{\rho' \omega'}_{\text{skew}} = \kappa_{\text{gm}} |S|^2 \overline{N}^2 \quad (92) $$

$$ \{C(E_k^t, E_a^t)\}_{\text{skew}} = \frac{g}{\rho_*} \overline{\rho' \omega'}_{\text{skew}} = -\kappa_{\text{gm}} |S|^2 \overline{N}^2 \quad (93) $$

Eq. (92) confirms that skew diffusion acts as a net sink of mean APE whose magnitude is proportional to the slope squared $|S|^2$ and mean squared buoyancy frequency $\overline{N}^2$ e.g., Griffies (2004).

### 4.4. Remarks on the eddy (turbulent) kinetic energy budget

Since skew diffusion associated with the meso-scale eddy parameterisation acts as a net source of eddy KE rather than of eddy APE, it is useful to conclude this section with some comments on the eddy KE budget and how the KE and APE cascades may help constrain mixing parameterisations. Under stationary and homogeneous conditions, the main sources of eddy KE are:

1. the downscale energy transfer associated with the mixing of momentum, primarily contributed by the vertical shear,

$$ C(E_k^m, E_k^t) \approx A_v \left( \frac{\partial \overline{\mathbf{u}}}{\partial z} \right)^2, \quad (94) $$

2. the downscale transfer associated with the meso-scale eddy parameterisation, which controls the conversion $C(E_a^m, E_a^t)$;

3. the loss of energy associated with diapycnal mixing, which controls $C(E_k^t, E_a^t)$;

4. the viscous dissipation rate, dissipating $E_k^t$ into background potential energy or ‘heat’.

Summing up all these contributions leads to the following balance:

$$ A_v \left( \frac{\partial \overline{\mathbf{u}}}{\partial z} \right)^2 + \kappa_{\text{gm}} |S|^2 \overline{N}^2 - K_{\rho} \overline{N}^2 \approx \varepsilon_k^t. \quad (95) $$

To show the dependence of this balance on the Richardson number, it is useful to divide this relation by $\overline{N}^2$, yielding

$$ A_v R_i^{-1} + \kappa_{\text{gm}} |S|^2 \approx K_{\rho} + \frac{\varepsilon_k^t}{\overline{N}^2} \quad (96) $$

where

$$ R_i = \overline{N}^2 \left( \frac{\partial \overline{\mathbf{u}}}{\partial z} \right)^{-2} \quad (97) $$

is the Richardson number. This relation shows that the different turbulent mixing parameters $K_v$, $\kappa_{\text{gm}}$, and $K_\rho$ are not independent from each other. Notably, by using the relation $K_\rho \approx \Gamma \varepsilon_k^t / \overline{N}^2$, the above relation can be written as

$$ A_v R_i^{-1} + \kappa_{\text{gm}} |S|^2 \approx \frac{(\Gamma + 1) \varepsilon_k^t}{\overline{N}^2} \quad (98) $$

and can be interpreted as a diagnostic energy balance potentially useful for predicting the turbulent kinetic energy dissipation rate, which is a crucial element of energetically consistent meso-scale eddy parameterisations such as in the GEOMETRIC framework (e.g., Marshall et al., 2012; Mak et al., 2018; Torres et al., 2023) or in Jansen et al. (2015)
and Jansen et al. (2019), for instance. In regions where the vertical shear can be predicted by thermal wind balance,

$$ \frac{\partial \bar{u}}{\partial z} \approx -\frac{g}{\rho_* f} \mathbf{k} \times \nabla \bar{\rho}, \qquad (99) $$

the TKE budget may also be written as

$$ \left( A_v \frac{\overline{N}^2}{f^2} + \kappa_{gm} \right) |S|^2 \overline{N}^2 \approx (\Gamma + 1) \epsilon_k^t. \qquad (100) $$

Eq. (100) shows that in this regime, the APE cascade converting mean APE into eddy APE behaves analogously to the KE cascade converting mean KE into eddy KE, similarly to Lindborg (2006)'s relations (3), provided that

$$ \kappa_{gm} \approx \Gamma A_v \frac{\overline{N}^2}{f^2} \qquad (101) $$

which appears to be closely related to the case discussed by Greatbatch and Lamb (1990). Determining to what extent the present results can inform further developments of such parameterisations is beyond the scope of this paper and will be addressed in a subsequent study.

## 5. Summary and discussion

In this study, we have revisited the mean/eddy decomposition theory for local Available Potential Energy (APE) density, focusing on its application to the characterisation of meso-scale eddies, small-scale irreversible mixing, and their interaction with the large-scale circulation. While irreversible mixing and eddy features occur at scales much smaller than the planetary-scale circulation, these processes play a crucial role in the global energy budget, necessitating careful consideration to ensure that forcing and dissipation terms are computed in a mutually consistent manner. Our approach addresses the consistency issues associated with the definition of localised forms of local APE density, which are often used to describe the energetics of these processes. Our analysis demonstrates that a consistent formulation should be based on partitioning the local APE density into mean and eddy components. The resulting eddy APE density measures forces against the local mean density profile, aligning with physical intuition. This formulation differs from heuristic localised forms of APE density, which ignore the dependence on the Lorenz reference state. Our exact eddy APE density retains some dependence on the Lorenz reference state, differing from the heuristic form by a factor of $\Lambda = \overline{N}^2/N_0^2$. Previous studies indicate that while fluid parcels are very close to their reference position in most of the ocean interior, where $\Lambda \approx 1$, this is generally not the case in the polar regions or the Gulf Stream area (e.g., Saenz et al., 2015; Tailleux, 2016, 2013a) where our results should be of most practical use.

Our approach is physically and conceptually simpler than that of Scotti and White (2014), while being mathematically equivalent, as discussed in detail by Tailleux (2024). The exact eddy APE budget is easier to interpret, as its net conversion with mean APE and eddy KE depends solely on the diapycnal component of the turbulent density flux, in contrast to the quasi-geostrophic (QG) approximated version. This formulation is valid for arbitrarily large departures from the Lorenz reference state, characterised by the parameters $\Lambda$ and the norm of the slope vector $|S|$. A fundamental difference between the two frameworks concerns the mean to eddy APE conversion, which can only be positive in the QG approximation, but which can potentially be also negative in the exact case. The exact eddy APE density behaves similarly to the heuristic localised form, despite being defined by different mathematical expressions. Both forms agree that the buoyancy involved in local turbulent stirring/mixing is the buoyancy anomaly defined relative to the mean density field, but they differ in how they define the displacement, $\zeta' = z - z_m$ versus $\zeta' = z_0(\rho) - z_0(\bar{\rho})$. This has important consequences for the study of meso-scale eddy APE and turbulent stratified mixing where the two differ significantly, notably

regarding the theory for mixing efficiency and the determination of turbulent vertical mixing diffusivity, as is expected to be the case in the polar regions, for instance.

Our approach represents a significant step towards a more realistic and accurate treatment of the APE budget, which is crucial for developing energetically consistent parameterisations and numerical ocean models. Importantly, such progress can only be achieved by the local theory of APE, highlighting the limitations of Winters et al. (1995)'s global APE framework. To fully assess the implications for estimating meso-scale eddy APE and the study of turbulent stratified mixing in the oceans, future work will involve reformulating the present framework using a Large Eddy Simulation (LES) spatial filter instead of Reynolds averaging and accounting for the nonlinearities of the equation of state for seawater. Preliminary results suggest that these nonlinearities can occasionally cause a loss of convexity for the local APE density, potentially corresponding to thermobarically unstable situations. The exact mean/eddy decomposition of the APE density, while simple, provides valuable insights into the turbulent APE cascade. Understanding how to rigorously apply this decomposition to finite-amplitude APE density is the first step towards a true multi-scale analysis, of which spectral analysis is one example. Further research is needed to develop a joint multi-scale analysis of APE density and kinetic energy (KE) that can provide more insights into how information about the Lorenz reference state affects different scales.

In conclusion, our study has clarified the relationship between heuristic and exact forms of localised eddy APE density and has demonstrated the importance of considering the dependence on the Lorenz reference state for a consistent treatment of the APE budget. This work lays the foundation for future research into the development of energetically consistent ocean mixing parameterisations and the multi-scale analysis of APE and KE in turbulent stratified flows.

### CRediT authorship contribution statement

Rémi Tailleux: Writing – review & editing, Writing – original draft, Methodology, Formal analysis, Conceptualization. Guillaume Roullet: Writing – review & editing, 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.

### References

Aiki, H., Zhai, X., Greatbatch, R.J., 2016. Energetics of the global ocean: the role of mesoscale eddies. In: Behera, S.K., Yamagata, T. (Eds.), Indo-Pacific Climate Variability and Predictability. World Scientific, pp. 109–134.

Andrews, D.G., 1981. A note on potential energy density in a stratified compressible fluid. J. Fluid Mech. 107, 227–236.

Arthur, R.S., Venayagamoorthy, S.K., Koseff, J.R., Fringer, O.B., 2017. How we compute n matters to estimates of mixing in stratified flows. J. Fluid Mech. 831 (R2).

Bleck, R., 1985. On the conversion between mean and eddy components of potential and kinetic energy in isentropic and isopycnic coordinates. Dyn. Atmos. Ocean. 9, 17–37.

Butler, E.D., Oliver, K.I.C., Gregory, J.M., Tailleux, R., 2013. The ocean’s gravitational potential energy budget in a coupled climate model. Geophys. Res. Lett. 40, 5417–5422.

Crocco, L., 1937. Eine neue stromfunktion für die erforschung der bewegung der gas mit rotation. Z. Angew. Math. Mech.

Dewar, W.K., McWilliams, J.C., 2019. On energy and turbulent mixing in the thermocline. J. Adv. Model. Earth Syst. 11, 578–596. http://dx.doi.org/10.1029/2018MS001502, Earth Systems, 11, 578–596.
Eckart, C., 1948. An analysis of stirring and mixing processes in incompressible fluids. J. Mar. Res. 7, 265-275.

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

Gargett, A.E., Holloway, G., 1984. Dissipation and diffusion by internal wave breaking. J. Mar. Res. 42, 15-27.

Gent, P., McWilliams, J.C., 1990. Isopycnal mixing in ocean circulation models. J. Phys. Oceanogr. 20, 150-155.

Gent, P., Willebrand, J., McDougall, T.J., McWilliams, J.C., 1995. Parameterizing eddy-induced tracer transports in ocean circulation models. J. Phys. Oceanogr. 25, 463-474.

Greatbatch, R.J., Lamb, K.G., 1990. On parameterizing vertical mixing of momentum in non-eddy resolving ocean models. J. Phys. Oceanogr. 20, 1634-1637.

Gregg, M.C., 2021. Ocean Mixing. Cambridge University Press.

Griffies, S.M., 1998. The Gent-McWilliams skew flux. J. Phys. Oceanogr. 28, 831-841. http://dx.doi.org/10.1175/1520-0485(1998)028%3C0831:TGMSF%3E2.0.CO;2

Griffies, S.M., 2004. Fundamentals of Ocean Models. Princeton University Press.

Griffies, S.M., Gnanadesikan, A., Pacanowski, R.C., Larichev, V.D., Dukowicz, J.K., Smith, R.D., 1998. Isoneutral diffusion in a z-coordinate ocean model. J. Phys. Oceanogr. 28, 805-830. http://dx.doi.org/10.1175/1520-0485(1998)028%3C0805: IDIAZC%3E2.0.CO;2

Holliday, D., McIntyre, M.E., 1981. On potential energy density in an incompressible, stratified fluid. J. Fluid Mech. 107, 221-225.

Jansen, M.F., Adcroft, A.J., Hallberg, R., Held, I.M., 2015. Parameterization of eddy fluxes based on a mesoscale energy budget. Ocean. Model. 92, 28-41. http://dx.doi.org/10.1016/j.ocemod.2015.05.007.

Jansen, M.F., Adcroft, A., Khani, S., Kong, H., 2019. Toward an energetically consistent, resolution aware parameterization of ocean mesoscale eddies. J. Adv. Model. Earth Syst. 11, 2844-2860. http://dx.doi.org/10.1029/2019MS001750.

Kucharski, F., 1997. On the concept of exergy and available potential energy. Q. J. R. Meteorol. Soc. 123, 2141-2156.

Kucharski, F., 2001. The interpretation of available potential energy as exergy applied to layers of a stratified atmosphere. Int. J. Exergy 1, 25-30.

Lindborg, E., 2006. The energy cascade in a strongly stratified fluid. J. Fluid Mech. 550, 207-242.

Lindborg, E., Brethouwer, G., 2008. Vertical dispersion by stratified turbulence. J. Fluid Mech. 614, 303-314.

Loose, N., Bachman, S., Grooms, I., Jansen, M., 2023. Diagnosing scale-dependent energy cycles in a high-resolution isopycnal ocean model. J. Phys. Oceanogr. 53, 157-176.

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

Luecke, C.A., Arbic, B.K., Bassette, S.L., Richman, J.G., Shriver, J.F., Alford, M.H., Smedstad, O.M., Timko, P.G., Trossman, D.S., Wallcraft, A.J., 2017. The global mesoscale eddy available potential energy field in models and observations. J. Geophys. Res. Ocean. 122, 9126-9143. http://dx.doi.org/10.1002/2017JC013126.

MacCready, P., Giddings, S.N., 2016. The mechanical energy budget of a regional ocean model. J. Phys. Oceanogr. 46, 2719-2733.

Mak, J., Maddison, J.R., Marshall, D.P., Munday, D.R., 2018. Implementation of a geometrically informed and energetically constrained mesoscale eddy parameterization in an ocean circulation model. J. Phys. Oceanogr. 48, 2363-2382.

Margules, M., 1903. Über die energie der stürme. Järnb. Zentralanst. Meteorol. Wien 40, 1-26.

Marquet, P., 1991. On the concept of exergy and available enthalpy: application to atmospheric energetics. Q. J. R. Met. Soc. 117, 449-475.

Marshall, D.P., Maddison, J.R., Berloff, P.S., 2012. A framework for parameterising eddy potential vorticity fluxes. J. Phys. Oceanogr. 42, 539-557.

Middleton, L., Fine, E.C., MacKinnon, J.A., Alford, M.H., Taylor, J.R., 2021. Estimating dissipation rates associated with double diffusion. Geophis. Res. Lett. 48, http://dx.doi.org/10.1029/2021GL092779.

Middleton, L., Taylor, J.R., 2020. A general criterion for the release of background potential energy through double diffusion. J. Fluid Mech. 893.

Munk, W.H., Wunsch, C., 1998. Abyssal recipes ii: Energetics of tidal and wind mixing. Deep-Sea Res. 45A 197, 7-2010.

Novak, L., Tailleux, R., 2018. On the local view of atmospheric available potential energy. J. Atmos. Sci. 75, 1891-1907. http://dx.doi.org/10.1175/JAS-D-17-0330.1.

Nycander, J., 2011. Energy conversion, mixing energy, and neutral surfaces with a nonlinear equation of state. J. Phys. Oceanogr. 41, 28-41.

Oakey, N.S., 1982. Determination of the rate of dissipation of turbulent energy from simultaneous temperature and velocity shear microstructure measurements. J. Phys. Oceanogr. 12, 256-271. http://dx.doi.org/10.1175/1520-0485(1982)012%3C0256: DOTROD%3E2.0.CO;2

Osborn, T.R., Cox, C.S., 1972. Oceanic fine structure. Geophys. Fluid Dyn. 3, 321-345.

Redi, M.H., 1982. Oceanic isopycnal mixing by coordinate rotation. J. Phys. Oceanogr. 12, 1154-1158.

Roullet, G., Capet, X., Maze, G., 2014. Global interior eddy available potential energy diagnosed from argo floats. Geophys. Res. Lett. 41, 1651-1656.

Roullet, G., Klein, P., 2009. Available potential energy diagnosis in a direct numerical simulation of rotating stratified turbulence. J. Fluid Mech. 624, 45-55.

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

Scotti, A., Beardsley, R., Butman, B., 2006. On the interpretation of energy and energy fluxes of nonlinear waves: an example from massachusetts bay. J. Fluid Mech. 561, 103-112.

Scotti, A., White, B., 2014. Diagnosing mixing in stratified turbulent flows with a locally defined available potential energy. J. Fluid Mech. 740, 114-135.

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

Smith, R.K., Montgomery, M.T., Zhu, H., 2005. Buoyancy in tropical cyclones and other rapidly rotating atmospheric vortices. Dyn. Atmos. Ocean. 40, 189-208.

Tailleux, R., 2009. On the energetics of turbulent stratified mixing, irreversible thermodynamics, boussinesq models, and the ocean heat engine controversy. J. Fluid Mech. 638, 339-382.

Tailleux, R., 2010. Entropy versus APE production: on the buoyancy power input in the oceans energy cycle. Geophys. Res. Lett. 37.

Tailleux, R., 2013a. Available potential energy and exergy in stratified fluids. Ann. Rev. Fluid Mech. 45, 35-58.

Tailleux, R., 2013b. 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., 2013c. Irreversible compressible work and available potential energy dissipation. Phys. Scr. T 155, 014033.

Tailleux, R., 2016. Generalized patched potential density and thermodynamic neutral density: Two new physically based quasi-neutral density variables for ocean water masses analyses and circulation studies. J. Phys. Oceanogr. 46, 3571-3584. http://dx.doi.org/10.1175/JPO-D-16-0072.1.

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., 2024. Negative available potential energy dissipation as the fundamental criterion for double diffusive instabilities. J. Fluid Mech. 994, http://dx.doi.org/10.1017/jfm.2024.647.

Tailleux, R., Dubos, T., 2024. A simple and transparent method for improving the energetics and thermodynamics of seawater approximations: Static energy asymptotics (SEA). Ocean. Model. 188, http://dx.doi.org/10.1016/j.ocemod.2024.102339.

Tailleux, R., Wolf, G., 2023. On the links between thermobaricity, available potential energy, neutral directions, buoyancy forces, potential vorticity, and lateral stirring in the ocean. http://dx.doi.org/10.48550/arXiv.2202.00456, unpublished, available on https://arxiv.org/abs/2202.00456.

Thorpe, A.J., Hoskins, B.J., Innocenti, V., 1989. The parcel method in a baroclinic atmosphere. J. Atmos. Sci. 46, 1274-1284.

Torres, R., Waldman, R., Mak, J., Séférrian, R., 2023. Global estimation of the eddy kinetic energy dissipation from a diagnostic energy balance. Geophis. Res. Lett. 50.

Vazsonyi, A., 1945. On rotational gas flows. Quart. Appl. Math. 3, 29-37.

von Storch, J.S., Eden, C., Fast, I., Haak, H., Hernandez-Deckers, D., Maier-Reimer, E., Marotzke, J., Stammer, D., 2012. An estimate of the lorenz energy cycle for the world ocean based on the 1/10° STORM/NCEP simulation. J. Phys. Oceanogr. 42, 2185-2205.

Winters, K.B., Lombard, P.N., Riley, J.J., d'Asaro, E.A., 1995. Available potential energy and mixing in density stratified fluids. J. Fluid Mech. 289, 115-128.

Wykes, M.S.D., Hughes, G., Dalziel, S.B., 2015. On the meaning of mixing efficiency for buoyancy-driven mixing in stratified turbulent flows. J. Fluid Mech. 781, 261-275.

Zemskova, V.E., White, B.L., Scotti, A., 2015. Available potential energy and the general circulation: partitioning wind, buoyancy forcing, and diapycnal mixing. J. Phys. Oceanogr. 45, 1510-1531.

Zilitinkevich, S.S., Elperin, T., Kleeorin, N., Rogachevskii, I., Esau, I., 2013. A hierarchy of energy- and flux (efb) turbulence closure models for stably-stratified geophysical flows. Bound.-Layer Meteorol. 146, 341-373.