On the Determination of the 3D Velocity Field in Terms of Conserved Variables in a Compressible Ocean

Abstract

Explicit expressions of the 3D velocity field in terms of the conserved quantities of ideal fluid thermocline theory, namely the Bernoulli function, density, and potential vorticity, are generalised in this paper to a compressible ocean with a realistic nonlinear equation of state. The most general such expression is the ‘inactive wind’ solution, an exact nonlinear solution of the inviscid compressible Navier–Stokes equation that satisfies the continuity equation as a consequence of Ertel’s potential vorticity theorem. However, due to the non-uniqueness of the choice of the Bernoulli function, such expressions are not unique and primarily differ in the magnitude of their vertical velocity component. Due to the thermobaric nonlinearity of the equation of state, the expression for the 3D velocity field of a compressible ocean is found to resemble its ideal fluid counterpart only if constructed using the available form of the Bernoulli function, the Bernoulli equivalent of Lorenz’s available potential energy (APE). APE theory also naturally defines a quasi-material, approximately neutral density variable known as the Lorenz reference density. This density variable, in turn, defines a potential vorticity variable that is minimally affected by thermobaric production, thus providing all the necessary tools for extending most results of ideal fluid thermocline theory to a compressible ocean.

1. Introduction

Although oceanic observations have increased dramatically in recent decades, especially since the development of the Argo floats program (https://argo.ucsd.edu, accessed on 7 March 2023), they still primarily constrain the density stratification rather than the velocity field. As a result, even in the most sophisticated ocean state estimate products [1,2,3,4], the realism of the simulated 3D oceanic velocity field $\mathbf{v}=(u,v,w)$ can rarely, if ever, be ascertained directly. Indeed, even when current meters are available, they usually contain high-frequency motions linked to tides and gravity waves that are not usually captured by numerical ocean models, thus complicating the comparison. In any case, comparisons, when possible, are only limited to the horizontal velocity field as the vertical velocity is generally too small to be directly measurable.

Since, for all practical purposes, the 3D oceanic velocity field is not a measurable parameter of the system, the development of accurate dynamical theories that predict how to infer it from available oceanic observations plays an especially important role in physical oceanography. It is, therefore, no surprise that the subject has a long history going back to the early days of the discipline, which, to date, has primarily relied on the use of the dynamic method [5]. As is well known, the dynamic method assumes geostrophic and hydrostatic balances and provides an explicit expression of the horizontal velocity field in terms of the vertical integrals of the horizontal density gradients relative to the velocity field at some reference level. Understanding how to specify the latter is one of the main challenges of the dynamic method and has given rise to many different approaches to tackle it. The simplest approach, which is most often encountered in the literature, assumes the existence of a level of no motion at which the geostrophic velocity vanishes (see [6] for a discussion of its limitations). That approach has since been superseded by the beta-spiral method of [7] or the inverse method(s) of [8,9], whose connection is discussed in [10]. Refs. [11,12,13,14] are examples of studies that further discuss these ideas. At some point, there was hope that the sea surface height (SSH) signal measured by satellite altimetry would be accurate enough to constrain the surface geostrophic velocity field. However, this requires a more accurate determination of the geoid than is presently available. Moreover, the fact that the SSH signal is contaminated by small-scale and high-frequency processes while also reflecting transient variations of the interior density field complicates filtering out only the relevant part needed by the dynamic method [15]. More recently, measurements of Argo float displacements have shown promise for constraining the geostrophic velocity at the floats’ parking depth, as recently demonstrated in [16,17]. Nevertheless, some limitations remain, as Argo float displacements are also a priori impacted by energetic small-scale and transient ageostrophic motions in addition to the geostrophic flow, which the authors sought to mitigate by averaging over the Argo period. Outside the Argo period, or when focusing on time snapshots, alternative/complementary approaches for constraining the unknown reference level are still needed.

Although the geostrophic approximation is generally accepted to represent an accurate theory for the horizontal velocity field $\mathbf{u}=(u,v)$, how best to predict the vertical velocity field w is, in contrast, much less understood and still actively debated. Physically, there are two main fundamental approaches to thinking about the vertical velocity field. The first approach, which, as far as we are aware, underlies the computation of the vertical velocity field in all existing numerical ocean models [18], is based on vertically integrating the continuity equation ${\nabla }_h·\mathbf{u}+w_z=0$. Doing so yields an expression for w controlled by the horizontal velocity divergence relative to the vertical velocity at some reference level $z_r$

$$w(x,y,z,t)=\underset{w_r}{\underbrace{w(x,y,z_r,t)}}-{\int }_{z_r}^z{\nabla }_h·\mathbf{u}\mathrm{d}{z}^{\prime },$$

(1)

where ${\nabla }_h$ is the horizontal nabla operator. As is well known, using the continuity equation in conjunction with geostrophy yields the celebrated linear Sverdrup balance $\beta v=f\partial w/\partial z$, which may be integrated to yield

$$w=w_r+\frac{\beta }{f}{\int }_{z_r}^{z}v\mathrm{d}{z}^{\prime },$$

(2)

(see [19]), where f is the Coriolis parameter and $\beta =df/dy$, with y denoting the latitude. The second approach is based on extracting w from the conservation equation $\partial C/\partial t+\mathbf{v}·\nabla C=\dot{C}$ of any conserved tracer C, viz.,

$$w={\left(\frac{\partial C}{\partial z}\right)}^{-1}\left(\dot{C}-\frac{\partial C}{\partial t}-\mathbf{u}·{\nabla }_{h}C\right),$$

(3)

where $\dot{C}$ represents the non-material sinks/sources of C. Both approaches come with important theoretical challenges. In the first approach, the difficulty stems from the fact that the horizontal divergence ${\nabla }_h·\mathbf{u}$ is often dominated by the ageostrophic component ${\mathbf{u}}_{ag}$ of the velocity, in which case (2) may become very inaccurate. In the second approach, the difficulty is that the time derivative $\partial C/\partial t$, and perhaps, to a lesser extent, the diabatic term $\dot{C}$, may both be important for predicting the large-scale vertical velocity field, yet are often hard (if at all possible) to constrain observationally. Existing approaches, therefore, can be regarded as attempts to mitigate these difficulties in some way. For instance, the well-known omega equation is obtained by artfully combining two prognostic equations (for vorticity and buoyancy, respectively) in a way that eliminates the time-derivative in each equation to formulate a diagnostic elliptic problem for the vertical velocity (see [20] for a recent implementation). Another approach of interest that shows how to infer the vertical velocity from individual moorings is that proposed in [21].

Physically, one key reason that makes the second approach potentially the most attractive is that if one ignores the difficulties associated with the time-dependent and diabatic terms $\partial C/\partial t$ and $\dot{C}$ in (3), knowledge of the geostrophic velocity field is generally sufficient enough to accurately estimate the horizontal advection term $-\mathbf{u}·{\nabla }_{h}C$. In contrast, knowledge of the ageostrophic velocity component ${\mathbf{u}}_{ag}$ is often needed for an accurate determination of the horizontal velocity divergence in the first approach, as mentioned above. As a result, there has been much interest in seeking to exploit the existence of quasi-material, conserved quantities to construct explicit expressions of the steady-state 3D velocity field. Thus, if $C_1$ and $C_2$ represent two independent conserved quantities, it is well known that in a steady state, the 3D velocity should lie at the intersection of two iso-surfaces of $C_1$ and $C_2$. Mathematically, this implies that the velocity field can be written in the form

$$\mathbf{v}=\lambda \nabla C_1\times \nabla C_2,$$

(4)

for some scalar field $\lambda$. Physically, the constraint determining $\lambda$ is that the horizontal component of (4) be either exactly or approximately geostrophic. The conserved quantities $C_1$ and $C_2$ then determine the vertical velocity according to $w=\mathbf{u}·{\mathbf{S}}_1$ or $w=\mathbf{u}·{\mathbf{S}}_2$ or both, where ${\mathbf{S}}_1$ and ${\mathbf{S}}_2$ are the horizontal ‘slope’ vectors

$${\mathbf{S}}_i=-{\left(\frac{\partial C_i}{\partial z}\right)}^{-1}\mathbf{u}·{\nabla }_h{C}_i,i=\{1,2\}.$$

(5)

Thus, depending on the conserved quantities considered, the relations $w=\mathbf{u}·{\mathbf{S}}_1$ and $w=\mathbf{u}·{\mathbf{S}}_2$ may or may not define the same vertical velocity. In practice, this means that the problem to be solved might be over-determined, and its resolution may require the use of least-squares methods.

Although (4) can, in principle, be implemented in practice using a variety of conserved tracers, three particular quantities deserve special attention owing to their theoretical and dynamical importance in oceanography, namely the Bernoulli function B, potential vorticity Q, and density $\rho$. Physically, this is because these quantities are the ones that underlie the ideal fluid thermocline equations [22] that have formed the basis for most ocean circulation theories, as well as those originally considered by Needler [23], who pioneered the use of (4). Unfortunately, such quantities have, so far, been unambiguously defined only for an ideal fluid, which significantly hinders our ability to evaluate the usefulness of (4) for the real compressible ocean. Although several studies have since extended Needler’s approach, e.g.,  [24,25,26,27], most discussions with a few exceptions, e.g., [28,29], have remained limited to the case of an ideal fluid in geostrophic balance. To make progress, this paper aims to generalise Needler’s approach to the case of a fully compressible ocean in order to systematically test its usefulness in the future (which is beyond the scope of this paper). Section 2 reviews the basic properties of the ideal fluid thermocline equations and introduces the concept of the ‘available’ Bernoulli function. Section 3 shows how to generalise ideal fluid thermocline theory to compressible seawater and shows that the results are related to [30]’s ‘inactive wind’ solution. Section 4 discusses the results and some perspectives.

2. Absolute Velocity Field Based On Ideal Fluid Thermocline Theory

2.1. Thermodynamic Form of Ideal Fluid Thermocline Equations

We begin by briefly reviewing the properties and structure of the ideal fluid thermocline equations [22] that underlie Needler’s determination of the absolute velocity field in terms of conserved properties [23]. Although most of the material is standard, an important new element is the introduction of the concept of the ‘available’ Bernoulli function, where the term ‘available’ is meant to parallel that of Lorenz’s ‘available potential energy’ [31]. The ideal fluid thermocline equations are written in the following form

$$f\mathbf{k}\times \rho \mathbf{v}+\nabla p=-\rho g\mathbf{k},$$

(6)

$$\nabla ·\mathbf{v}=0,$$

(7)

$$\mathbf{v}·\nabla \rho =0,$$

(8)

and describe the geostrophic and hydrostatic momentum balances, continuity, and conservation of density, respectively, where p is the pressure, $\rho$ is the density, g is the acceleration of gravity, and $\mathbf{k}$ is the unit vector pointing upward. Although the geostrophic and hydrostatic balances are expected to meaningfully describe the balanced part of the flow even under transient evolution for a sufficiently small Rossby number, the steady form of the density equation in (8) is, in contrast, much more difficult to justify rigorously, especially without a clarification of what the density variable $\rho$ is supposed to represent. Presumably, any justification of (8) must consider some form of temporal averaging that does not introduce eddy-correlation terms, such as thickness-weighted averaging [32]. As the issue is quite complex and cannot fully be addressed without having first clarified the exact nature of $\rho$, its full treatment is deferred to a subsequent study.

The key step for linking the velocity field to the conserved quantities of the system is to rewrite the momentum in (6) in its thermodynamic form (also known as the Crocco–Vazsonyi form [33,34]),

$$f\mathbf{k}\times \rho \mathbf{v}+{\rho }_{\star }\nabla B_h^{ideal}={\rho }_{\star }{\mathbf{P}}_h^{ideal},{\mathbf{P}}_h^{ideal}=\frac{gz\nabla \rho }{{\rho }_{\star }},$$

(9)

where $B_h^{ideal}=(p+\rho gz)/{\rho }_{\star }$ is the standard Bernoulli function and ${\rho }_{\star }$ is a constant Boussinesq reference density.

2.2. Available and Background Bernoulli Functions

One of the main ideas of this paper is that Lorenz’s APE theory holds the key to understanding how to generalise the key ingredients of ideal fluid thermocline theory, that is, the Bernoulli function, potential vorticity, and density, to a compressible ocean. The key ingredients needed here are the reference pressure and density profiles, denoted by $p_0\left(z\right)$ and ${\rho }_0\left(z\right)$, respectively. These profiles characterise the Lorenz reference state of minimum potential energy that can be obtained through an adiabatic and isohaline re-arrangement of mass. Additionally, we require the reference position of a fluid parcel, $z_r$, which is defined as the solution of the level of neutral buoyancy (LNB) equation $\rho ={\rho }_0\left(z_r\right)$. This definition yields $z_r=z_r\left(\rho \right)$ as a function of density only [35,36,37]. Physically, the APE theory is useful for distinguishing between the potential energy that can be converted into kinetic energy in a reversible way (the APE) and the dynamically inert part of potential energy (the BPE) that cannot. Because the Bernoulli function $B_h^{ideal}=(p+\rho gz)/{\rho }_{\star }$ is thermodynamic in nature, the same idea that only a fraction of it is available for reversible conversions with kinetic energy must apply. This motivates us to define its dynamically inert part in a Lagrangian sense as its value in a Lorenz reference state, viz.,

$$B_r^{ideal}=\frac{p_0\left(z_r\right)+{\rho }_0\left(z_r\right)g{z}_r}{{\rho }_{\star }},$$

(10)

and the ‘available’ part of the Bernoulli function as 

$$\begin{array}{cc}\hfill {\rho }_{\star }{B}_a^{ideal}={\rho }_{\star }(B_h^{ideal}-B_r^{ideal})=& p-p_0\left(z\right)+p_0\left(z\right)-p_0\left(z_r\right)+\rho g(z-z_r)\hfill \\ \hfill =& p-p_0\left(z\right)+\underset{{\rho }_{\star }{E}_a}{\underbrace{{\int }_{z_r}^{z}g(\rho -{\rho }_0\left(z^{\prime }\right))\mathrm{d}{z}^{\prime }}}\hfill \end{array}$$

(11)

using the fact that, by definition, $d{p}_0/dz\left(z\right)=-g{\rho }_0\left(z\right)$ and ${\rho }_0\left(z_r\right)=\rho$, where the underbraced quantity $E_a$ can be recognised as the positive definite APE density that was originally introduced in [38,39] and later extended to multi-component Boussinesq and stratified fluids in [35,37]. Physically, the APE density can be written as the work against buoyancy forces needed to bring a fluid parcel from its level of neutral buoyancy $z_r$ (hence, satisfying $b(S,\theta ,z_r)=0\left)\right)$ to its actual position, viz.,

$$E_a=-{\int }_{z_r}^{z}b(S,\theta ,z^{\prime })\mathrm{d}{z}^{\prime },$$

(12)

and the only difference between the expressions for a Boussinesq and a general compressible fluid is in the expressions for the buoyancy b. For a Boussinesq fluid, $b=-g(\rho -{\rho }_0\left(z\right))/{\rho }_{\star }$, whereas for compressible seawater, which is discussed in the next section, $b=-g[1-{\rho }_0\left(z\right)\upsilon (S,\theta ,p_0\left(z\right))]$. Note that since $\rho ={\rho }_0\left(z_r\right)$, ${\rho }_{\star }{E}_a$ may also be rewritten as

$${\rho }_{\star }{E}_a=-g{\int }_{z_r}^z{\int }_{z_r}^{z^{\prime }}\frac{d{\rho }_0}{dz}\left(z^{″}\right)\mathrm{d}{z}^{″}\mathrm{d}{z}^{\prime },$$

(13)

so that for a small departure from the reference position,

$$E_a\approx -\frac{g}{{\rho }_{\star }}\frac{d{\rho }_0}{dz}\left(z_r\right)\frac{{(z-z_r)}^2}{2}=\frac{N_r^2{(z-z_r)}^2}{2}$$

(14)

which a reader unfamiliar with APE may still recognise. It also makes the positive definite character of $E_a$ clearer if needed.

Removing $\nabla B_r$ from both sides of (9) using the easily verified result that $\nabla B_r=g{z}_r\nabla \rho$, yields the following available thermodynamic form of momentum balance

$$f\mathbf{k}\times \rho \mathbf{v}+{\rho }_{\star }\nabla B_a^{ideal}={\rho }_{\star }{\mathbf{P}}_a^{ideal},{\rho }_{\star }{\mathbf{P}}_a^{ideal}=g(z-z_r)\nabla \rho .$$

(15)

A key point to note here is that in both (9) and (15), the two P-vectors ${\mathbf{P}}_h^{ideal}$ and ${\mathbf{P}}_a^{ideal}$ are parallel to the gradient of density $\nabla \rho$, and, therefore, perpendicular to the isopycnal surfaces $\rho =\mathrm{constant}$. We see in Section 3 that this property is lost in a compressible ocean. It is of interest to note that ${\mathbf{P}}_h^{ideal}$ is the ideal fluid counterpart of the P-vector previously identified in [40].

2.3. Bernoulli and Potential Vorticity (PV) Theorems

The Bernoulli theorem [41] and Ertel’s PV conservation theorem [42,43] play key roles in this paper. The former is trivially obtained by taking the inner product of (9) and (15) with $\mathbf{v}$, accounting for (8), which immediately yields $\mathbf{v}·\nabla B_h=0$ and $\mathbf{v}·\nabla B_a=0$. The proof of the potential vorticity (PV) conservation theorem is somewhat more involved. To obtain it, first, divide (6) by $\rho$ and take the curl, thus leading to

$$\mathbf{k}(\mathbf{v}·\nabla f)-f\frac{\partial \mathbf{v}}{\partial z}=\frac{\nabla \rho \times \nabla p}{{\rho }^2}.$$

(16)

Next, take the inner product of (16) with $\nabla \rho$, which then yields

$$\frac{\partial \rho }{\partial z}\mathbf{v}·\nabla f-f\nabla \rho ·\frac{\partial \mathbf{v}}{\partial z}=\mathbf{v}·\nabla \left(f\frac{\partial \rho }{\partial z}\right)=0,$$

(17)

which establishes the material conservation of the potential vorticity (PV) $Q=(f/{\rho }_{\star })\partial \rho /\partial z$ as expected.

2.4. Determination of the Absolute Velocity Field in Terms of Conserved Quantities

Taking the cross product of $\nabla \rho$ with (9) or (15) and using the result that $\nabla \rho \times (f\mathbf{k}\times \mathbf{v})=(\mathbf{v}·\nabla \rho )f\mathbf{k}-(f\nabla \rho ·\mathbf{k})\mathbf{v}=-f{\rho }_z\mathbf{v}$ (which follows from the vector algebra relation $\mathbf{A}\times (\mathbf{B}\times \mathbf{C})=(\mathbf{A}·\mathbf{C})\mathbf{B}-(\mathbf{A}·\mathbf{B})\mathbf{C}$ and conservation of density $\mathbf{v}·\nabla \rho =0$, with ${\rho }_z$ being shorthand for $\partial \rho /\partial z$) leads, after some manipulation, to the following explicit expression of $\mathbf{v}$ that was previously obtained in [23],

$$\rho \mathbf{v}=\frac{\nabla \rho \times \nabla B}{Q},$$

(18)

which explicitly relies on density and the Bernoulli function being conserved following fluid parcels, regardless of which form of the Bernoulli function is used (B refers indifferently to $B_h$ or $B_a$).

To show that (18) naturally satisfies the continuity equation as a consequence of Ertel’s PV conservation theorem, simply take its divergence, which yields

$$\nabla ·\left(\rho \mathbf{v}\right)=-\frac{(\nabla \rho \times \nabla B)·\nabla Q}{Q^2}=-\frac{\rho \mathbf{v}·\nabla Q}{Q}=0,$$

(19)

QED. As a result, (18) represents an exact steady solution of the ideal fluid Equations (6)–(8) that, in addition to the horizontal velocity, predicts the vertical velocity

$$\rho w=\frac{\mathbf{k}·({\nabla }_h\rho \times {\nabla }_{h}B)}{Q}.$$

(20)

Importantly, Needler’s formula (22) is insensitive to the particular choice of Bernoulli function—B or $B_a$—used to estimate it due to $B_r=B_r\left(\rho \right)$ being a function of $\rho$ only.

2.5. Bernoulli Method

Physically, one key reason that Needler’s formula (18) is appealing is that the idea that the steady or time-averaged 3D velocity field should lie at the intersection of the iso-surfaces of two conserved quantities is a priori valid beyond the geostrophic approximation. Its other key advantage is that it naturally satisfies the continuity equation as a consequence of Ertel’s potential vorticity conservation theorem, which also generalises beyond the geostrophic approximation, as discussed in the next section. In the context of the dynamic method, however, Needler’s formula does not, in itself, solve the problem of the unknown reference level. Indeed, because ${\nabla }_{h}B=\nabla p_h+gz{\nabla }_h\rho$, it is easily verified that the horizontal velocity it predicts is the standard geostrophic balance, whereas the vertical velocity predicted in (20) is predicted more easily from the density equation $\mathbf{v}·\rho =0$, i.e.,

$$w=-{\left(\frac{\partial \rho }{\partial z}\right)}^{-1}\mathbf{u}·{\nabla }_h\rho =\mathbf{u}·{\mathbf{S}}_{\rho },$$

(21)

where ${\mathbf{S}}_{\rho }$ is the slope vector associated with the isopycnal surfaces. In other words, Needler’s formula does not solve the problem of the unknown reference level because the Bernoulli function depends on the same unknown constant of integration as the pressure field.

To circumvent the difficulty and make (18) useful, it is necessary to invoke the result that if $\rho$, B, and Q are all conserved along fluid parcel trajectories, they must be functionally related. Thus, to the extent that such a result holds and this functional relationship may be written in the form $B=G(\rho ,Q)$, it is possible to rewrite (18) as 

$$\rho \mathbf{v}=\frac{\partial G}{\partial Q}\frac{\nabla \rho \times \nabla Q}{Q},$$

(22)

[23]. Now, the key advantage of  (22) over (18) is that the term $(\nabla \rho \times \nabla Q)/Q$ can, in principle, be empirically estimated from climatological fields without having to solve the problem of the unknown reference level since the latter affects neither $\rho$ nor Q. In other words, (22) provides a determination of the absolute velocity field up to the multiplicative constant $\partial G/\partial Q$. In this approach, the problem of the unknown reference level, therefore, transforms into the problem of how best to evaluate the functional relationship $G(\rho ,Q)$ and the partial derivative $\partial G/\partial Q$. For related discussions of these ideas, the reader is referred to, e.g., [13,24,25,26,27,28,44,45]. However, it seems fair to say that the current empirical evidence for a well-defined relationship $B=B(\rho ,Q)$ is inconclusive at best. Current approaches, however, have relied on using the conventional form of the Bernoulli function $B=B_h$ rather than its available form $B=B_a$. It will be of interest in future studies to test whether using $B_a$, as well as the density variable discussed in the next section, can lead to a better-defined relationship $B_a=B_a(\rho ,Q)$.

3. Generalisation to Compressible Seawater

3.1. Governing Equations for Compressible Seawater

We now turn to a realistic nonlinear ocean described by the compressible Navier–Stokes equations, treating seawater as a two-constituent stratified fluid:

$$\frac{D\mathbf{v}}{Dt}+\mathbf{2}\mathrm{\Omega }\times \mathbf{v}+\frac{1}{\rho }\nabla p=-\nabla \mathrm{\Phi }+\mathbf{F},$$

(23)

$$\frac{\partial \rho }{\partial t}+\nabla ·\left(\rho \mathbf{v}\right)=0,$$

(24)

$$\frac{D\eta }{Dt}=\dot{\eta },\frac{DS}{Dt}=\dot{S},$$

(25)

$$\upsilon =\frac{1}{\rho }=\upsilon (\eta ,S,p),$$

(26)

where $\eta$ is the specific entropy, S is the salinity, $\upsilon =1/\rho$ is the specific volume, $\mathbf{F}$ is a friction force, and $\mathrm{\Phi }=gz$ is the geopotential. Moreover, $\dot{\eta }$ and $\dot{S}$ denote the diabatic sources/sinks of heat and salt due to molecular diffusive fluxes (radiation can be included if needed). As in the previous section, the first step is to rewrite the momentum balance (23) in its thermodynamic or Crocco–Vazsonyi form

$$\frac{\partial \mathbf{v}}{\partial t}+{\omega }_a\times \mathbf{v}+\nabla B_h={\mathbf{P}}_h+\mathbf{F},$$

(27)

where ${\omega }_a=\xi +2\mathrm{\Omega }$ is the absolute vorticity and $\xi =\nabla \times \mathbf{v}$ is the relative vorticity. The results were obtained using the well-known relation $(\mathbf{v}·\nabla )\mathbf{v}=\nabla ({\mathbf{v}}^2/2)+\xi \times \mathbf{v}$, as well as the total differential for the specific enthalpy $\mathrm{d}h=T\mathrm{d}\eta +\mu \mathrm{d}S+\upsilon \mathrm{d}p$, where T is the in situ temperature and $\mu$ is the relative chemical potential. For a compressible ocean, the conventional form of the Bernoulli function is

$$B_h=\frac{{\mathbf{v}}^2}{2}+h(\eta ,S,p)+\mathrm{\Phi }\left(z\right),$$

(28)

whereas the P-vector ${\mathbf{P}}_h$ is the vector that retains the thermohaline gradient of h, viz.,

$${\mathbf{P}}_h=\frac{\partial h}{\partial \eta }\nabla \eta +\frac{\partial h}{\partial S}\nabla S=T\nabla \eta +\mu \nabla S.$$

(29)

As in the previous section, we introduce the available Bernoulli function as $B_a=B_h-B_r$, that is, as the difference between the conventional form of the Bernoulli function $B_h$ and the background Bernoulli function

$$B_r=h(\eta ,S,p_0\left(z_r\right))+\mathrm{\Phi }\left(z_r\right)$$

(30)

where, as before, $z_r=z_r(\eta ,S)$ is the reference height of a fluid parcel in Lorenz’s reference state of minimum potential energy, which, in practice, may be computed using the computationally efficient algorithm of Saenz et al. [36]. Physically, $z_r$ is defined as before as a solution of the level of neutral buoyancy equation, which, for two-component compressible seawater, takes the form

$$\upsilon (\eta ,S,p_0\left(z_r\right))={\upsilon }_0\left(z_r\right),$$

(31)

where ${\upsilon }_0\left(z\right)=1/{\rho }_0\left(z\right)$. The available Bernoulli function may thus be written in the form

$$\begin{array}{cc}\hfill B_a=& \frac{{\mathbf{v}}^2}{2}+h(\eta ,S,p)-h(\eta ,S,p_r)+g(z-z_r)\hfill \\ \hfill =& \frac{{\mathbf{v}}^2}{2}+\mathrm{\Pi }+\frac{p-p_0\left(z\right)}{\rho }\hfill \end{array}$$

(32)

where $\mathrm{\Pi }={\mathrm{\Pi }}_1+{\mathrm{\Pi }}_2$ is the potential energy density defined in [37], and ${\mathrm{\Pi }}_1$ and ${\mathrm{\Pi }}_2$ are the subcomponents

$${\mathrm{\Pi }}_1=h(\eta ,S,p)-h(\eta ,S,p_0\left(z\right))+\frac{p_0\left(z\right)-p}{\rho },$$

(33)

$${\mathrm{\Pi }}_2=h(\eta ,S,p_0\left(z\right))-h(\eta ,S,p_0\left(z_r\right))+g(z-z_r)$$

(34)

As explained and demonstrated in [37], both ${\mathrm{\Pi }}_1$ and ${\mathrm{\Pi }}_2$ are positive definite. ${\mathrm{\Pi }}_1$ may be interpreted as the available compressible energy (ACE), which represents the expansion/contraction work required to compress/expand from the reference pressure $p_0\left(z\right)$ to the actual pressure. As for ${\mathrm{\Pi }}_2$, it is the APE density and represents the work against the buoyancy forces required to move the fluid parcel from its reference position $z_r$ at pressure $p_0\left(z_r\right)=p_r$ to its actual position z at pressure $p_0\left(z\right)$. As before, the available thermodynamic form of the momentum balance is obtained by removing $\nabla B_r=T_r\nabla \eta +{\mu }_r\nabla S$ from both sides of (27), which leads to

$$\frac{\partial \mathbf{v}}{\partial t}+{\omega }_a\times \mathbf{v}+\nabla B_a={\mathbf{P}}_a+\mathbf{F}$$

(35)

where the P-vector ${\mathbf{P}}_a$ takes the form

$${\mathbf{P}}_a=(T-T_r)\nabla \eta +(\mu -{\mu }_r)\nabla S.$$

(36)

Note that the suffix ‘r’ denotes the thermodynamic quantities estimated at the reference pressure $p_r=p_0\left(z_r\right)$.

3.2. Comparison of Ideal and Compressible Forms of Bernoulli Functions and P-Vectors

In an ideal fluid, Needler’s formula (18) defines the 3D velocity field as being perpendicular to the density surfaces $\rho =\mathrm{constant}$, regardless of the particular form of $(B,\mathbf{P})$ used due to ${\mathbf{P}}_h^{ideal}$ and ${\mathbf{P}}_a^{ideal}$ both being parallel to $\nabla \rho$ in that case. The 3D velocity field is also perpendicular to the iso-surfaces of the Bernoulli function in both cases but, generally, the iso-surfaces of $B_h^{ideal}$ and $B_a^{ideal}$ should appear distinct from each other. Whether the available form $(B_a^{ideal},{\mathbf{P}}_a^{ideal})$ is superior to the conventional form $(B_h^{ideal},{\mathbf{P}}_h^{ideal})$ cannot be determined from theoretical considerations alone; however, it is possible that the assumed empirical functional relationship $B=B(\rho ,Q)$ underlying the Bernoulli method may be more accurately achieved in climatological observations of the density field by using one of the forms of $(B,\mathbf{P})$.

The situation is different for compressible seawater, however, as ${\mathbf{P}}_a$ and ${\mathbf{P}}_h$ now, generally, define different directions, as is clear from Equations (29) and (36). Moreover, it is easily verified that both P-vectors have a non-zero helicity $H_p=\mathbf{P}·(\nabla \times \mathbf{P})$ because of thermobaricity and are, therefore, both non-integrable. Mathematically, this means that it is not possible to identify a well-defined seawater variable whose iso-surfaces are everywhere perpendicular to $\mathbf{P}$. From a theoretical perspective, this is both interesting and important, as it suggests that the superior form of $(B,\mathbf{P})$ may be the one for which the compressible and ideal expressions resemble each other the most.

To explore this idea,

Table 1. Comparison of the ideal and compressible forms of the Bernoulli function and P-vectors. Note that the Bernoulli functions for a compressible fluid do not include the kinetic energy term.

Quantity	Ideal	Compressible
$B_h$	$\frac{p+\rho gz}{{\rho }_{\star }}$	$h(\eta ,S,p)+gz$
$B_a$	$E_a+\frac{p-p_0\left(z\right)}{{\rho }_{\star }}$	${\mathrm{\Pi }}_1+{\mathrm{\Pi }}_2+\frac{p-p_0\left(z\right)}{\rho }$
${\mathbf{P}}_h$	$\frac{gz\nabla \rho }{{\rho }_{\star }}$	$T\nabla \eta +\mu \nabla S$
${\mathbf{P}}_a$	$\frac{g(z-z_r)\nabla \rho }{{\rho }_{\star }}$	$(T-T_r)\nabla \eta +(\mu -{\mu }_r)\nabla S$

summarises the forms of the Bernoulli function and P-vectors for the ideal and compressible forms established in the previous sections, whereas Figure 1 illustrates the four different types of Bernoulli functions along the ${30}^{\circ }$ W meridional section in the Atlantic Ocean. We used the WOCE climatological dataset [46] and computed the hydrostatic pressure, assuming a level of no motion at 1500 $\mathrm{m}$. The Lorenz reference density and pressure profiles, as well as the reference depths, were computed as in [47]. The contribution from the kinetic energy was ignored. The figure clearly shows that the Bernoulli function depends sensitively on the approach considered, as well on the choice of the arbitrary constants entering the definition of the specific enthalpy, as defined by TEOS-10 (www.teos-10.org, accessed on 7 March 2023). In the figure, panel (a) shows that the ideal standard Bernoulli function is dominated by its depth variations, with values increasing in height from −700 $\mathrm{J}.{\mathrm{kg}}^{-1}$ at depth to close to 0 at the surface. In contrast, the compressible standard Bernoulli function (panel (b)), which depends on the TEOS-10 definition of the specific enthalpy, is dominated by thermal variations $\propto c_{p0}\theta$ near the surface that result in much bigger values overall, which also increase in height from about close to zero at depth to about 60,000 $\mathrm{J}.{\mathrm{kg}}^{-1}$ near the surface. The ideal and compressible forms of the available Bernoulli function are depicted in panels (c) and (d), respectively. Although our theory predicts that the latter should be close to each other, this is not the case in practice because our prediction in (11) assumes an incompressible $\rho$, which is not true of in situ density. Panel (d) shows that the compressible available Bernoulli function $B_a$ is dominated by horizontal variations. This seems to be an optimal behaviour for plotting it on isopycnal surfaces, which will be discussed in detail in a subsequent study. Note that the values exhibited by $B_a$, which range approximately from 0 to 20 $\mathrm{J}.{\mathrm{kg}}^{-1}$, are considerably smaller than that of all other Bernoulli functions, confirming the idea that $B_a$ is the only one not affected by irrelevant dynamical information, which we hope can be exploited in the future.

3.2.1. Conventional Bernoulli Function

By contrasting the ideal and compressible forms of the conventional Bernoulli function,

$$\frac{p+\rho gz}{{\rho }_{\star }}\leftrightarrow h(\eta ,S,p)+gz$$

(37)

it is clear that the two expressions are not easily related. The difference in behaviour between $B_h$ and $B_h^{ideal}$ can be further evidenced by contrasting, for instance, their vertical derivatives, leading to

$$\frac{\partial }{\partial z}{B}_h^{ideal}=\frac{\partial }{\partial z}\left(\frac{p+\rho gz}{{\rho }_{\star }}\right)=\frac{gz}{{\rho }_{\star }}\frac{\partial \rho }{\partial z},$$

(38)

$$\frac{\partial B_h}{\partial z}=\frac{\partial }{\partial z}(h+gz)=T\frac{\partial \eta }{\partial z}+\mu \frac{\partial S}{\partial z}.$$

(39)

In order to link $\partial B_h/\partial z$ to the vertical derivative of some density variable, one needs to make use of the Maxwell relations attached to the total differential of the enthalpy $\mathrm{d}h=T\mathrm{d}\eta +\mu \mathrm{d}S+\upsilon \mathrm{d}p$, which is obtained by stating the equality of the cross derivatives, viz.,

$$\frac{{\partial }^{2}h}{\partial \eta \partial p}=\frac{\partial T}{\partial p}=\frac{\partial \upsilon }{\partial \eta }\frac{{\partial }^{2}h}{\partial S\partial p}=\frac{\partial \mu }{\partial p}=\frac{\partial \upsilon }{\partial S}.$$

(40)

To proceed, let us now introduce a reference pressure $p_{\star }=p_0\left(z_{\star }\right)$ that is envisioned as being not too different from p, defining the reference values $T_{\star }=T(\eta ,S,p_{\star })$ and ${\mu }_{\star }=\mu (\eta ,S,p_{\star })$. If we also use use the approximation $p\approx p_0\left(z\right)$, (39) can be rewritten as 

$$\begin{array}{cc}\hfill \frac{\partial B_h}{\partial z}=& {\int }_{p_{\star }}^p\frac{\partial \upsilon }{\partial \eta }\mathrm{d}{p}^{\prime }\frac{\partial \eta }{\partial z}+{\int }_{p_{\star }}^p\frac{\partial \upsilon }{\partial S}\mathrm{d}{p}^{\prime }\frac{\partial S}{\partial z}+T_{\star }\frac{\partial \eta }{\partial z}+{\mu }_{\star }\frac{\partial S}{\partial z}\hfill \\ \hfill \approx & (p-p_{\star })\left(\overline{{\upsilon }_{\eta }}\frac{\partial \eta }{\partial z}+\overline{{\upsilon }_S}\frac{\partial S}{\partial z}\right)+T_{\star }\frac{\partial \eta }{\partial z}+{\mu }_{\star }\frac{\partial S}{\partial z}\hfill \\ \hfill \approx & \frac{g(z-z_{\star })}{{\rho }_{\star }}\frac{\partial {\rho }_{lr}}{\partial z}+T_{\star }\frac{\partial \eta }{\partial z}+{\mu }_{\star }\frac{\partial S}{\partial z}.\hfill \end{array}$$

(41)

where ${\rho }_{lr}$ is the potential density referenced to $p_{\star }$ (or more accurately to $(p_{\star }+p)/2$), where the result assumes that $p_{\star }$ is close enough to p that the overbar value can be approximated by their values at $p_{\star }$. Although the first term of (41) succeeds in showing dependence on the vertical derivative of some density variable that makes it directly comparable to $d{B}_h^{ideal}/dz$, the second term is extraneous. Using the same ideas, it is also easily established that the same difficulties exist in relating ${\mathbf{P}}_h$ to ${\mathbf{P}}_h^{ideal}$.

3.2.2. Available Bernoulli Function

In contrast, it is immediately apparent that the ideal and compressible available Bernoulli functions

$$E_a+\frac{p-p_0\left(z\right)}{{\rho }_{\star }}\leftrightarrow {\mathrm{\Pi }}_1+{\mathrm{\Pi }}_2+\frac{p-p_0\left(z\right)}{\rho }$$

(42)

are directly comparable, as $B_a^{ideal}$ can be obtained from $B_a$ by neglecting ${\mathrm{\Pi }}_1$ and replacing $\rho$ with ${\rho }_{\star }$ in $(p-p_0\left(z\right))/\rho$, as Tailleux [37] showed that ${\mathrm{\Pi }}_2$ closely resembled $E_a$. Likewise, it is possible to show that the APE-based P-vectors resemble each other

$$\frac{g(z-z_r)\nabla \rho }{{\rho }_{\star }}\leftrightarrow (T-T_r)\nabla \eta +(\mu -{\mu }_r)\nabla S$$

(43)

To see the connection, one may use the Maxwell relations (40) to establish that

$$\begin{array}{cc}\hfill {\mathbf{P}}_a=& (T-T_r)\nabla \eta +(\mu -{\mu }_r)\nabla S\hfill \\ \hfill =& {\int }_{p_r}^p\frac{\partial \upsilon }{\partial \eta }(\eta ,S,p^{\prime })\mathrm{d}{p}^{\prime }\nabla \eta +{\int }_{p_r}^p\frac{\partial \upsilon }{\partial S}(\eta ,S,p^{\prime })\mathrm{d}{p}^{\prime }\nabla S\hfill \\ \hfill =& (p-p_r)\left(\overline{{\upsilon }_{\eta }}\nabla \eta +\overline{{\upsilon }_S}\nabla S\right)\hfill \\ \hfill \approx & \frac{g(z-z_r)}{\overline{\rho }}\left(\overline{{\rho }_{\eta }}\nabla \eta +\overline{{\rho }_S}\nabla S\right)\hfill \end{array}$$

(44)

using the fact that ${\upsilon }_{\eta }=-{\rho }_{\eta }/{\rho }^2$, ${\upsilon }_S=-{\rho }_S/{\rho }^2$, and $p-p_r\approx -\overline{\rho }g(z-z_r)$, where ${\upsilon }_{\eta }=\partial \upsilon /\partial \eta$, ${\upsilon }_S=\partial \upsilon /\partial S$, etc. This shows that ${\mathbf{P}}_a$ is intermediate between the locally-referenced gradients referenced to the local and reference pressure, respectively, which means that it is close to the standard neutral vector considered in [48]. In other words, ${\mathbf{P}}_a$ is approximately parallel to the gradient of the potential density referenced to the mid pressure $(p+p_r)/2$. For a more extensive discussion of this link, see [49].

The above considerations clearly establish that the results of ideal fluid thermocline theory can only be generalised to compressible seawater if the available form $(B_a,{\mathbf{P}}_a)$ rather than the conventional form $(B_h,{\mathbf{P}}_h)$ is used.

3.3. Inactive Wind Solutions

To discuss how to extend Needler’s formula to compressible seawater, we consider the steady and inviscid momentum equations written in their thermodynamic, or Crocco–Vazsonyi, form

$${\omega }_a\times \mathbf{v}+\nabla B_{\ell }={\mathbf{P}}_{\ell },$$

(45)

(from Equations (27) and (35)), where the subscript $\ell =a,h$ indicates whether the conventional or available form of B and $\mathbf{P}$ is used. To simplify the notations, this subscript is dropped in the following and re-introduced only when needed. For a small Rossby number, as pertains to the large-scale motions of interest here, the relative vorticity and kinetic energy only affect ${\omega }_a$ and B at second order. Equation (45) is then an under-determined linear system for $\mathbf{v}$ that can only be solved if the following solvability condition is satisfied

$${\omega }_a·\mathbf{P}={\omega }_a·\nabla B,$$

(46)

(obtained by taking the inner product of (45) with ${\omega }_a$), so that $\mathbf{v}$ in (45) is determined only up to an arbitrary vector $\chi {\omega }_a$ parallel to ${\omega }_a$, with $\chi$ a scalar field. The $\chi =0$ solution of (45) that is perpendicular to both $\nabla B$ and $\mathbf{P}$ and is the counterpart of Needler’s formula (18), is easily verified to be

$$\mathbf{v}\stackrel{\mathrm{def}}{=}{\mathbf{v}}_{ia}=\frac{\mathbf{P}\times \nabla B}{\rho Q_b},$$

(47)

whereas the counterpart of $Q=f{\rho }_z$ is

$$Q_b=\frac{{\omega }_a·\nabla B}{\rho }=\frac{{\omega }_a·\mathbf{P}}{\rho },$$

(48)

which is the equality following from (46). In the context of a dry atmosphere, [30] derived an expression similar to (47) and referred to ${\mathbf{v}}_{ia}$ as the ‘inactive wind’.

In order to examine the consistency of the inactive wind solution (47) with mass conservation, let us take the divergence of $\rho {\mathbf{v}}_{ia}$, which yields

$$\begin{array}{cc}\hfill \nabla ·\left(\rho {\mathbf{v}}_{ia}\right)=& -\frac{(\mathbf{P}\times \nabla B)}{Q_b^2}·\nabla Q_b+\frac{(\nabla \times \mathbf{P})·\nabla B}{Q_b}\hfill \\ \hfill =& -\frac{\rho {\mathbf{v}}_{ia}·\nabla Q_b}{Q_b}+\frac{(\nabla \times \mathbf{P})·\nabla B}{Q_b}.\hfill \end{array}$$

(49)

Similarly as for one of Needler’s formulas, ${\mathbf{v}}_{ia}$ is shown here to satisfy the continuity equation as a consequence of Ertel’s PV conservation theorem applied to the PV constructed from the Bernoulli function $Q={\omega }_a·\nabla B/\rho$. To show this, let us recall that in its most general form, [50]’s theorem (see [43] for an English translation) establishes that for any scalar $\lambda$, the PV variable $Q_{\lambda }={\omega }_a·\nabla \lambda /\rho$ can be shown to satisfy the conservation law

$$\frac{D{Q}_{\lambda }}{Dt}=\frac{{\omega }_a}{\rho }·\nabla \left(\frac{D\lambda }{Dt}\right)+\frac{1}{{\rho }^3}\nabla \lambda ·(\nabla \rho \times \nabla p)+\frac{\nabla \lambda ·\nabla \times \mathbf{F}}{\rho }$$

(50)

(e.g., see Equation 4.95 in [51]). To make the link with (49), it is important to understand the different equivalent forms that the baroclinic production term (the term proportional to $\nabla \rho \times \nabla p$ in (50)) may assume. As seen previously, from the definitions of ${\mathbf{P}}_h$ and ${\mathbf{P}}_a$, we have the following equivalence relations

$$\frac{1}{\rho }\nabla p+\nabla \mathrm{\Phi }=\nabla (h+\mathrm{\Phi })-{\mathbf{P}}_h=\nabla (h+\mathrm{\Phi }-B_r)-{\mathbf{P}}_a.$$

(51)

Taking the curl yields the following equivalent expressions

$$\frac{1}{{\rho }^2}\nabla \rho \times \nabla p=\frac{1}{{\rho }^2}\mathbf{N}\times \nabla p=\nabla \times {\mathbf{P}}_h=\nabla \times {\mathbf{P}}_a.$$

(52)

which are proportional to the baroclinic production term, where $\mathbf{N}={\rho }_S\nabla S+{\rho }_{\theta }\nabla \theta ={\rho }_{\eta }\nabla \eta +{\rho }_S\nabla S$ is the so-called N-neutral vector entering [48]’s definition of (approximately) neutral surfaces (ANS). Now, if we use $B=\lambda$ in the inviscid case $\mathbf{F}=0$ and invoke the Bernoulli theorem, $DB/Dt=0$, (50) predicts that

$$\frac{D{Q}_b}{Dt}={\mathbf{v}}_{ia}·\nabla Q_b=\frac{1}{\rho }\nabla B·(\nabla \times \mathbf{P})$$

(53)

Comparing this with (49) shows that

$$\nabla ·\left(\rho {\mathbf{v}}_{ia}\right)=0,$$

(54)

as expected, thus confirming that ${\mathbf{v}}_{ia}$ satisfies the continuity equation regardless of which form of $(B,\mathbf{P})$ is used.

3.4. Uniqueness of the Inactive Wind Solution

The two inactive wind solutions introduced above may be written explicitly as

$${\mathbf{v}}_{ia}^h=\frac{{\mathbf{P}}_h\times \nabla B_h}{{\omega }_a·\nabla B_h},{\mathbf{v}}_{ia}^a=\frac{{\mathbf{P}}_a\times \nabla B_a}{{\omega }_a·\nabla B_a}.$$

(55)

Since these two solutions both satisfy the continuity equation, a question that naturally arises is whether they define the same velocity field despite being constructed from different fields. To examine this, let us recall that, by definition, ${\mathbf{P}}_a={\mathbf{P}}_h-\nabla B_r$, whereas $B_a=B_h-B_r$. It follows that

$$\begin{array}{cc}\hfill {\mathbf{P}}_a\times \nabla B_a=& ({\mathbf{P}}_h-\nabla B_r)\times \nabla (B_h-B_r)\hfill \\ \hfill =& {\mathbf{P}}_h\times \nabla B_h+(\nabla B_h-{\mathbf{P}}_h)\times \nabla B_r\hfill \\ \hfill =& {\mathbf{P}}_h\times \nabla B_h+\left(\nabla \frac{{\mathbf{v}}^2}{2}+\frac{1}{\rho }{\nabla }_{h}p\right)\times \nabla B_r,\hfill \end{array}$$

(56)

which suggests that ${\mathbf{v}}_{ia}^h$ and ${\mathbf{v}}_{ia}^a$ define two different vector fields. However, because they both represent an exact solution of (45), it follows that their difference must be proportional to the null-space solution ${\omega }_a$. In other words, there must exist some scalar $\delta \lambda$ such that

$${\mathbf{v}}_{ia}^h-{\mathbf{v}}_{ia}^a=\delta \lambda {\omega }_a.$$

(57)

For a small Rossby number, ${\omega }_a\approx f\mathbf{k}$ (where f is the Coriolis parameter), which implies that ${\mathbf{v}}_{ia}^h$ and ${\mathbf{v}}_{ia}^a$ must primarily differ by their vertical velocity component. Physically, this is plausible because in this case, the horizontal components of ${\mathbf{v}}_{ia}^h$ and ${\mathbf{v}}_{ia}^a$ must both be approximately geostrophic, i.e., ${\mathbf{u}}_{ia}^h\approx {\mathbf{u}}_{ia}^a\approx {\mathbf{u}}_g$, whereas their vertical components must satisfy $w_{ia}^h={\mathbf{u}}_g·{\mathbf{S}}_h$ and $w_{ia}^a={\mathbf{u}}_g·{\mathbf{S}}_a$, respectively, where ${\mathbf{S}}_h$ and ${\mathbf{S}}_a$ are the horizontal slope vectors defined by ${\mathbf{P}}_h$ and ${\mathbf{P}}_a$, respectively. As a result,

$$w_{ia}^h-w_{ia}^a={\mathbf{u}}_g·({\mathbf{S}}_h-{\mathbf{S}}_a),$$

(58)

which confirms that the two vertical velocities $w_{ia}^h$ and $w_{ia}^a$ might differ if the slopes defined by the two different P-vectors differ, as is generally the case. It is of interest to ascertain whether the vertical velocity component of ${\mathbf{v}}_{ia}^a$ is a better predictor of the actual vertical velocity field than that of ${\mathbf{v}}_{ia}^h$, which we plan on investigating in a subsequent study.

3.5. Reformulation in Terms of Quantities Independent of Pressure

Let ${\gamma }^T={\gamma }^T(S,\theta )$ denote an approximately neutral quasi-material density variable and $Q_{\gamma T}={\omega }_a·\nabla {\gamma }^T/\rho$ denote the PV variable constructed from it. As discussed in [36,47,49,52], APE theory naturally includes a generalised form of the potential density ${\rho }^{LZ}(S,\theta )=\rho (S,\theta ,p_r)$ that is naturally very accurately neutral outside the Southern Ocean (even more so than Jackett and McDougall [53]’s empirical neutral density variable ${\gamma }^n)$ while also being mathematically and physically well defined. This motivated [47,52] to define thermodynamic neutral density ${\gamma }^T$ as 

$${\gamma }^T=\rho (S,\theta ,p_r)-f_n\left(p_r\right)$$

(59)

where $f_n\left(p_r\right)$ is a polynomial pressure correction empirically fitted to make ${\gamma }^T$ look as much like ${\gamma }^n$. Because it tends to be more accurately neutral than ${\gamma }^n$ outside the Southern Ocean, ${\gamma }^T$ is the variable that is currently the least affected by the thermobaric production term.

If one accepts ${\gamma }^T$ and $Q_{\gamma T}$ as the most sensible generalisations of the concepts of density and PV to compressible seawater, then one may proceed similarly to [23] and assume that if B, ${\gamma }^T$, and $Q_{\gamma T}$ are all approximately conserved along fluid parcel trajectories, a functional relationship should exist between B, ${\gamma }^T$, and $Q_{\gamma T}$, say $B=G({\gamma }^T,Q_{\gamma T})$, as in ideal fluid thermocline theory. Then, the seawater counterpart of  Needler’s second formula (22) becomes

$${\mathbf{v}}_{ia}=\frac{\partial G}{\partial {\gamma }^T}\frac{\mathbf{P}\times \nabla {\gamma }^T}{\rho Q_b}+\frac{\partial G}{\partial Q_{\gamma T}}\frac{\mathbf{P}\times \nabla Q_{\gamma T}}{\rho Q_b}.$$

(60)

In comparison to Needler’s formula, (60) possesses the extra and undesirable thermobaricity-induced term proportional to $\mathbf{P}\times \nabla {\gamma }^T$. This term can only be neglected if ${\mathbf{P}}_a$ rather than ${\mathbf{P}}_h$ is used, as the angle between ${\mathbf{P}}_h$ and $\nabla {\gamma }^T$ is not generally small enough. If so,

$${\mathbf{v}}_{ia}\approx \frac{\partial G}{\partial Q_{\gamma T}}\frac{{\mathbf{P}}_a\times \nabla Q_{\gamma T}}{\rho Q_b},$$

(61)

which is more directly comparable to (22), stressing again the fundamental importance of APE theory to extend the results of ideal fluid thermocline theory to compressible seawater.

4. Discussion

The idea that steady fluid parcel trajectories lie at the intersection of the iso-surfaces of conserved quantities is arguably one of the most promising avenues of research for progressing the theory of the 3D oceanic velocity field. This is because this idea is a priori as equally valid for an ideal fluid for which it was originally developed by  [23] as it is for a fully compressible ocean with a realistic nonlinear equation of state, as recently initiated by [29]. In this paper, we made significant progress towards generalising this idea to compressible seawater, encapsulated into three main new results.

Our first main result is that Needler’s formula (18) can be interpreted as a linear approximation of a much more general nonlinear and exact solution of the compressible NSE called the ’inactive wind’ solution. This solution was previously derived by [30] in the context of a dry atmosphere and was extended here to two-component compressible seawater. Like Needler’s formula, the inactive wind solution satisfies the continuity equation as a consequence of Ertel’s PV conservation theorem [42] (that is itself related to the generalised Bernoulli theorem of [41]). Like Needler’s formula, the inactive wind solution is perpendicular to the gradient of the Bernoulli function but unlike Needler’s formula, it is perpendicular to a vector $\mathbf{P}$ rather than to the gradient of density $\nabla \rho$. Physically, the inactive wind solution is most easily obtained by rewriting the momentum equations in their Crocco–Vazsonyi, or thermodynamic, form, which appears to be the most illuminating form for relating the 3D velocity field to the conserved quantities of the fluid.

Our second main result is that both Needler’s formula and the inactive wind solution are sensitive to how the Bernoulli function is defined, as it is always possible to redefine the latter by subtracting an arbitrary quasi-material function of $\rho$ (for a simple fluid) or S and $\theta$ (for compressible seawater). We find that only if the ‘available’ form of the Bernoulli function is used is it possible to meaningfully relate the inactive wind solution to Needler’s formula. Physically, the available Bernoulli function is defined as the difference between the conventional Bernoulli function and its background reference value in Lorenz’s state of minimum potential energy entering Lorenz’s APE theory [31]. Indeed, only in this case is the vector $\mathbf{P}$ parallel to an approximately neutral density variable, namely the Lorenz reference density (LRD), which is discussed at length in [47,49,52].

Our third main result is that inactive wind solutions defined for different forms of the Bernoulli function and vector $\mathbf{P}$ do not necessarily define the same 3D velocity field, even if each represents an exact solution of the compressible NSE satisfying the continuity equation. Mathematically, this is because the nonlinear balance equation from which the inactive wind solution is determined is degenerate. For a small Rossby number, this is equivalent to saying that different inactive wind solutions differ primarily in their vertical velocity components. More generally, this result means that although in principle, it is possible to construct a 3D velocity field as $\lambda \nabla C_1\times \nabla C_2$ in terms of any arbitrary conserved quantities $C_1$ and $C_2$ for some $\lambda$, this does not necessarily imply that all constructions define the same velocity field, which does not appear to have received much attention so far. Physically, this is important because it provides the means, at least in principle, to test the usefulness of different constructions by comparing the 3D velocity field that each expression predicts with the 3D velocity field from any dynamically consistent ocean state estimate, as we plan on pursuing in a subsequent study.

The present results are important because we believe that they can pave the way for a more rigorous and general theory of the oceanic 3D velocity field valid for a realistic compressible ocean, as we hope to further demonstrate through concrete applications in future studies.