On the determination of the 3D velocity field in terms of conserved variables in a compressible ocean

TL;DR

This paper generalizes the explicit expressions for the 3D velocity field in a compressible ocean using conserved variables, incorporating thermobaric effects and extending ideal fluid thermocline theory.

Contribution

It introduces a generalized 'inactive wind' solution for compressible oceans and connects it with Lorenz APE theory and potential vorticity, enhancing understanding of ocean dynamics.

Findings

Expressions for 3D velocity fields are non-unique due to Bernoulli function choices.

The velocity field resembles ideal fluid solutions when constructed with Lorenz APE theory.

Lorenz reference density provides a minimally affected potential vorticity variable.

Abstract

Explicit expressions of the 3D velocity field in terms of the conserved quantities of ideal fluid thermocline theory, namely Bernoulli function, density, and potential vorticity, are generalised here 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 that satisfies the continuity equation as a consequence of Ertel's potential vorticity theorem. Such expressions are shown to be non-unique due to the non-uniqueness of the choice of Bernoulli function and in general approximately differ by the magnitude of their vertical velocity component. Due to the thermobaric nonlinearity of the equation of state, the expression of the 3D velocity field for a compressible ocean is found to resemble its ideal fluid counterpart only if constructed in terms of the available form of Bernoulli function as per Lorenz theory of available potential energy (APE). APE theory also naturally defines a quasi-material approximately neutral density variable called Lorenz reference density, which in turn defines a potential vorticity variable minimally affected by thermobaric production, thus providing all necessary tools for extending most results of ideal fluid thermocline theory to a compressible ocean.

Tables

Table 1. Table 1: Comparison of the ideal and compressible forms of Bernoulli function and P-vectors. Note that the Bernoulli functions for a compressible fluid don’t 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 }}$	${\Pi }_1+{\Pi }_2+\frac{p-p_0\left(z\right)}{\rho }$
${𝐏}_h$	$\frac{gz\nabla \rho }{{\rho }_{\star }}$	$T\nabla \eta +\mu \nabla S$
${𝐏}_a$	$\frac{g\left(z-z_r\right)\nabla \rho }{{\rho }_{\star }}$	$\left(T-T_r\right)\nabla \eta +\left(\mu -{\mu }_r\right)\nabla S$

Equations

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

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

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

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

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

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

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

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

$$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 }},$$

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

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

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

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

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

$$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 .$$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

$$B_r=h(\eta ,S,p_0(z_r))+\Phi (z_r)$$

Full text

\pubvolume

1 \issuenum1 \articlenumber0

\datereceived \daterevised \dateaccepted \datepublished \hreflinkhttps://doi.org/ \TitleOn the determination of the 3D velocity field in terms of conserved variables in a compressible ocean \TitleCitationOn the determination of the 3D velocity field in terms of conserved variables in a compressible ocean \AuthorRémi Tailleux 1\orcidA \AuthorNamesRémi Tailleux \AuthorCitationTailleux, R.

\corresCorrespondence: R.G.J.Tailleux@reading.ac.uk

1 Introduction

Although oceanic observations have increased dramatically over the past decades, especially since the development of the Argo floats program (https://argo.ucsd.edu), these still primarily constrain the density stratification rather than the velocity field. As a result, even in the most sophisticated ocean state estimate products Lee et al. (2009); Forget et al. (2015); Masuda and Osafune (2021); Wunsch and Ferrari (2019), the realism of the simulated 3D oceanic velocity field ${\bf v}=(u,v,w)$ can rarely if ever be ascertained directly. Indeed, even when current meters are available, these usually contain high-frequency motions linked to tides and gravity waves 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 predicting 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 Fomin (1964). 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 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 and most often encountered one in the literature assumes the existence of a level of no motion at which the geostrophic velocity vanishes, see Olbers and Willebrand (1984) for a discussion of its limitations. Since then, such an approach has been superseded by the beta-spiral method of Schott and Stommel (1978) or the inverse method(s) of Wunsch (1977, 1978), whose connection is discussed in Davis (1978). Behringer (1979); Behringer and Stommel (1980); Killworth (1986); Bigg (1985) are examples of studies further discussing 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 Park (2004). More recently, measurements of Argo floats displacements have shown promise for constraining the geostrophic velocity at the floats’ parking depth, as recently demonstrated by Colin de Verdière and Ollitrault (2016); Colin de Verdière et al. (2019). Nevertheless, some limitations remain as Argo floats 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.

While the geostrophic approximation is generally accepted to represent an accurate theory for the horizontal velocity field ${\bf 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 for 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 Griffies (2004), is based on vertically integrating the continuity equation $\nabla_{h}\cdot{\bf 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}$

[TABLE]

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

[TABLE]

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

[TABLE]

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\cdot{\bf u}$ is often dominated by the ageostrophic component ${\bf u}_{ag}$ of the velocity, in which case (2) may become very inaccurate. In the second approach, the difficult stems from 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 often hard (if at all possible) to constrain observationally. Existing approaches, therefore, can be regarded as attempts to mitigate such 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 eliminate the time-derivative in each equation to formulate a diagnostic elliptic problem for the vertical velocity, see Nardelli (2020) for a recent implementation. Another approach of interest showing how to infer the vertical velocity from individual moorings is that proposed by Sevellec et al. (2015).

Physically, one key reason that makes the second approach potentially the most attractive one is because 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 for accurately estimating the horizontal advection term $-{\bf u}\cdot\nabla C$. In contrast, knowledge of the ageostrophic velocity component ${\bf 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, then 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 may be written in the form

[TABLE]

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

[TABLE]

Thus, depending on the conserved quantities considered, the relations $w={\bf u}\cdot{\bf S}_{1}$ and $w={\bf u}\cdot{\bf 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 that its resolution may need the use of least-squares.

While (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 underlying the ideal fluid thermocline equtions Welander that have formed the basis for most ocean circulation theories, as well as those originally considered by Needler Needler (1985), who pioneered the use of (4). Unfortunately, such quantities have been so far unambiguously defined only for an ideal fluid, which significantly hinder our ability to evaluate the usefulness of (4) for the real compressible ocean. Although several studies have extended Needler’s approach since, e.g., Chu (1995, 2000); Kurgansky et al. (2002); Kurgansky (2021), most discussions with a few exceptions, e.g., McDougall (1995); Ochoa et al. (2020) have remained limited to the case of an ideal fluid in geostrophic balance. To make progress, this paper aims to show how to generalise Needler’s approach to the case of a fully compressible ocean, in order to be able to test its usefulness systematically in the future (which is beyond the scope of this paper). Section 2 reviews the basic properties of the ideal fluid thermocline equations and introduce the concept of ’available’ Bernoulli function. Section 3 shows how to generalise ideal fluid thermocline theory to compressible seawater, and shows that the results is related to Gassmann (2014) ‘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 Welander that underlie Needler’s determination of the absolute velocity field in terms of conserved properties Needler (1985). Although most of the material is standard, an important new element is the introduction of the concept of ‘available’ Bernoulli function, where the term ‘available’ is meant to parallel that of Lorenz ‘available potential energy’ Lorenz (1955). The ideal fluid thermocline equations are here written in the following form

[TABLE]

[TABLE]

[TABLE]

and respectively describe the geostrophic and hydrostatic momentum balance, continuity, and conservation of density respectively, where ${\bf v}=(u,v,w)$ is the three-dimensional velocity field, $p$ is pressure, $\rho$ is density, $g$ is the acceleration of gravity, $f$ is the Coriolis parameter, and ${\bf k}$ is the unit vector pointing upward. While the geostrophic and hydrostatic balances are expected to meaningfully describe the balanced part of the flow even under transient evolution for sufficiently small Rossby number, the steady form of the density equation (8) is in contrast much more difficult to justify rigorously, especially without clarity as to what density variable $\rho$ is supposed to represent. Presumably, any justification of (8) must require consideration of some form of temporal averaging that does not introduce eddy-correlation terms, such as thickness-weighted averaging Young (2012). 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 main key step for linking the velocity field to the conserved quantities of the system is to rewrite the momentum (6) in its thermodynamic form (also known as the Crocco-Vazsonyi form (Crocco, 1937; Vazsonyi, 1945)),

[TABLE]

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

A key idea of this paper is that Lorenz 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 $p_{0}(z)$ and $\rho_{0}(z)$ characterising Lorenz reference state of minimun potential energy obtainable by means of an adiabatic an isohaline re-arrangement of mass, and the reference position of a fluid parcel $z_{r}$ defined as a solution of the level of neutral buoyancy (LNB) equation $\rho=\rho_{0}(z_{r})$, which defines $z_{r}=z_{r}(\rho)$ as a function of density only Tailleux (2013); Saenz et al. (2015); Tailleux (2018). Physically, APE theory is useful for providing a physical mean to tease apart the part of the potential energy that is available for reversible conversions with kinetic energy (the APE) from the dynamically inert part (the BPE) that is not. 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 Lorenz reference state, viz.,

[TABLE]

and the ‘available’ part of the Bernoulli function as

[TABLE]

using the fact that $p_{0}^{\prime}(z)=-g\rho_{0}(z)$ and $\rho_{0}(z_{r})=\rho$ by definition, where the underbraced quantity $E_{a}$ can be recognised as the positive definite APE density originally introduced by Holliday and McIntyre (1981) and Andrews (1981) and later extended to multi-component Boussinesq and stratified fluids by Tailleux (2013, 2018). Physically, the APE density may 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))$ to its actual position, viz.,

[TABLE]

the only difference between the expressions for a Boussinesq and general compressible fluid being in the expressions for the buoyancy $b$. For a Boussinesq fluid, $b=-g(\rho-\rho_{0}(z))/\rho_{\star}$, while for compressible seawater discussed in next section $b=-g[1-\rho_{0}(z)\upsilon(S,\theta,p_{0}(z))]$. Note that since $\rho=\rho_{0}(z_{r})$, $\rho_{\star}E_{a}$ may also be rewritten as

[TABLE]

so that for small departure from the reference position,

[TABLE]

which the reader unfamiliar with APE may still recognise, and which makes the positive definite character of $E_{a}$ clearer if needed.

If we remove $\nabla B_{r}$ from both sides of (9), using the easily verified result that $\nabla B_{r}=gz_{r}\nabla\rho$, yields the following available thermodynamic form of momentum balance

[TABLE]

A key point to note here is that in both (9) and (15), the two P-vectors ${\bf P}_{h}^{ideal}$ and ${\bf P}_{a}^{ideal}$ are proportional to the gradient of density $\nabla\rho$, and therefore perpendicular to the isopycnal surfaces $\rho={\rm constant}$. We shall see in Section 3 that this property is lost in a compressible ocean. It is of interest to note that ${\bf P}_{h}^{ideal}$ is the ideal fluid counterpart of the P-vector previously identified by Nycander (2011).

2.3 Bernoulli and Potential vorticity (PV) theorems

The Bernoulli theorem Schär (1993) and Ertel’s PV conservation theorem Muller (1995); Schubert et al. (2004) play a key role in this paper. As regards to the former, it is trivially obtained by taking the inner product of (9) and (15) by ${\bf v}$, accounting for (8), which immediately yields ${\bf v}\cdot\nabla B_{h}=0$ and ${\bf v}\cdot\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

[TABLE]

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

[TABLE]

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{\bf k}\times{\bf v})=({\bf v}\cdot\nabla\rho)f{\bf k}-(f\nabla\rho\cdot{\bf k}){\bf v}=-f\rho_{z}{\bf v}$ (which follows from the vector algebra relation ${\bf A}\times({\bf B}\times{\bf C})=({\bf A}\cdot{\bf C}){\bf B}-({\bf A}\cdot{\bf B}){\bf C}$ and conservation of density ${\bf v}\cdot\nabla\rho=0$, with $\rho_{z}$ being shorhand for $\partial\rho/\partial z$) leads after some manipulation to the following explicit expression of ${\bf v}$ previously obtained by Needler (1985),

[TABLE]

which explicitly relies on density and the Bernoulli function being conserved following fluid parcels, regardless of which form of 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

[TABLE]

QED. As a result, (18) represents an exact steady solution of the ideal fluid equations (6-8) that also predicts the vertical velocity

[TABLE]

not just the horizontal velocity, where $\nabla_{h}$ denotes the horizontal gradient. 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}(\rho)$ being a function of $\rho$ only.

2.5 Bernoulli method

Physically, one key reason why Needler’s formula (18) is appealing is because 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 potential vorticity conservation theorem, which also generalises beyond the geostrophic approximation, as discussed in 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 that it predicts is the standard geostrophic balance, while the vertical velocity predicted by (20) is actually the one predicted more easily from the density equation ${\bf v}\cdot\nabla\rho=0$, i.e.,

[TABLE]

where ${\bf 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 turn (18) into something concretely 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 that this functional relationship may be written in the form $B=G(\rho,Q)$, it become possible to rewrite (18) as

[TABLE]

Needler (1985). 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 neither affects $\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., Killworth (1979, 1980, 1986); Kurgansky et al. (2002); Kurgansky (2021); Chu (1995, 2000); McDougall (1995). 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 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 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:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\eta$ is the specific entropy, $S$ is salinity, $\upsilon=1/\rho$ is the specific volume, ${\bf F}$ is a friction force, and $\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

[TABLE]

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

[TABLE]

while the P-vector ${\bf P}_{h}$ is the vector that retains the thermohaline gradient of $h$, viz.,

[TABLE]

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 Bernoulli function $B_{h}$ and the background Bernoulli function

[TABLE]

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

[TABLE]

where $\upsilon_{0}(z)=1/\rho_{0}(z)$. The available Bernoulli function may thus be written in the form

[TABLE]

where $\Pi=\Pi_{1}+\Pi_{2}$ is the potential energy density defined in Tailleux (2018), with $\Pi_{1}$ and $\Pi_{2}$ the subcomponents

[TABLE]

[TABLE]

As explained and demonstrated in Tailleux (2018), both $\Pi_{1}$ and $\Pi_{2}$ are positive definite. $\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}(z)$ to the actual pressure. As to $\Pi_{2}$, it is the APE density and represents the work against buoyancy forces required to move the fluid parcel from its reference position $z_{r}$ at pressure $p_{0}(z_{r})=p_{r}$ to its actual position $z$ at pressure $p_{0}(z)$. As before, the available thermodynamic form of momentum balance is obtained by removing $\nabla B_{r}=T_{r}\nabla\eta+\mu_{r}\nabla S$ on both sides of (27), which leads to

[TABLE]

where the P-vector ${\bf P}_{a}$ takes the form

[TABLE]

Note that the suffix ’r’ denote thermodynamic quantities estimated at the reference pressure $p_{r}=p_{0}(z_{r})$.

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={\rm constant}$ regardless of which particular form of $(B,{\bf P})$ is used, due to ${\bf P}_{h}^{ideal}$ and ${\bf P}_{a}^{ideal}$ being both 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 the iso-surfaces of $B_{h}^{ideal}$ and $B_{a}^{ideal}$ should in general appear distinct from each other. Whether the available form $(B_{a}^{ideal},{\bf P}_{a}^{ideal})$ is superior over the conventional form $(B_{h}^{ideal},{\bf 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 might be more accurately achieved in climatological observations of the density field for one of $(B,{\bf P})$ forms.

The situation is different for compressible seawater, however, as ${\bf P}_{a}$ and ${\bf P}_{h}$ now define different directions in general, as is clear from their expressions (29) and (36). Moreover, it is easily verified that both P-vectors have a non zero helicity $H_{p}={\bf P}\cdot(\nabla\times{\bf 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 ${\bf P}$. From a theoretical viewpoint, this is both interesting and important, as it suggests that the superior form of $(B,{\bf P})$ may be the one for which the compressible and ideal expressions resemble the most to each other.

To explore this idea, table 1 summarises the forms of Bernoulli function and P-vectors for an ideal and compressible forms established in the previous sections, while Fig. 1 illustrates the 4 different types of Bernoulli functions along the $30^{\circ}W$ meridional section in the Atlantic Ocean. We used the WOCE climatological dataset Gouretski and Koltermann (2004) and computed the hydrostatic pressure assuming a level of no motion at $1500\,{\rm m}$. The Lorenz reference density and pressure profiles, as well as the computation of the reference depths, were computed as in Tailleux (2021). The contribution from the kinetic energy was ignored. This figure clearly shows that the Bernoulli function depends sensitively on the approach considered as well as on the choice of the arbitrary constants entering the definition of the specific enthalpy as defined by TEOS-10 (www.teos-10.org). Panel (a) shows the ideal standard Bernoulli function to be dominated by its depth variations, with values increasing with height from $-700\,{\rm J.kg^{-1}}$ at depth to close to [math] at the surface. In contrast, the compressible standard Bernoulli function (Panel (b)), which depends on the TEOS-10 definition of specific enthalpy, is dominated by thermal variations $\propto c_{p0}\theta$ near the surface that results in much bigger values overall, which also increase with height from about close to zero at depth to about $60,000\,{\rm J.kg^{-1}}$ near the surface. The ideal and compressible forms of Bernoulli function are depicted in Panels (c) and (d) respectively. Although our theory predicts that the ideal and compressible available Bernoulli functions should be close to each other, this is not the case in practice because our prediction (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, which seems an optimal behaviour for plotting it on isopycnal surfaces as will be discussed in details in a subsequent study. Note that the range of values exhibited by $B_{a}$, comprised between 0 and $20\,{\rm J.kg^{-1}}$ approximately, is 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.

Conventional Bernoulli function

Contrasting the ideal and compressible forms of the conventional Bernoulli function,

[TABLE]

makes it 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 their vertical derivative for instance, leading to

[TABLE]

[TABLE]

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 enthalpy ${\rm d}h=T{\rm d}\eta+\mu{\rm d}S+\upsilon{\rm d}p$, obtained by stating the equality of the cross derivatives, viz.,

[TABLE]

To proceed, let us now introduce some reference pressure $p_{\star}=p_{0}(z_{\star})$ 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}(z)$, (39) may be rewritten as

[TABLE]

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

Available Bernoulli function

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

[TABLE]

are directly comparable, as $B_{a}^{ideal}$ can be obtained from $B_{a}$ by neglecting $\Pi_{1}$ and replacing $\rho$ by $\rho_{\star}$ in $(p-p_{0}(z))/\rho$, as Tailleux (2018) showed $\Pi_{2}$ to closely resemble $E_{a}$. Likewise, it is possible to show that the APE-based P-vectors resemble each other

[TABLE]

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

[TABLE]

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 result shows that ${\bf 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 by McDougall (1987). In other words, ${\bf P}_{a}$ is approximately proportional to the gradient of the potential density referenced to the mid-pressure $(p+p_{r})/2$. For more extensive discussion of this link, see Tailleux and Wolf (2023).

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},{\bf P}_{a})$ rather than the conventional form $(B_{h},{\bf 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

[TABLE]

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

[TABLE]

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

[TABLE]

while the counterpart of $Q=f\rho_{z}$ is

[TABLE]

the equality following from (46). In the context of a dry atmosphere, Gassmann (2014) derived an expression similar to (47) and referred to ${\bf v}_{ia}$ as the inactive wind.

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

[TABLE]

As for one of Needler’s formula, it can be shown that ${\bf v}_{ia}$ satisfies the continuity equation as a consequence of Ertel’s PV conservation theorem applied to the PV constructed from the Bernoulli function $Q=\boldsymbol{\omega}_{a}\cdot\nabla B/\rho$. To show this, let us recall that in its most general form, Ertel (1942)’s theorem (see Schubert et al. (2004) for an English translation) establishes that for any scalar $\lambda$, the PV variable $Q_{\lambda}=\boldsymbol{\omega}_{a}\cdot\nabla\lambda/\rho$ can be shown to satisfy the conservation law

[TABLE]

e.g., see Eq. (4.95) of Vallis (2006). 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 ${\bf P}_{h}$ and ${\bf P}_{a}$, we have the following equivalence relations

[TABLE]

Taking the curl yields the following equivalent expressions

[TABLE]

which are proportional to the baroclinic production term, where ${\bf 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 McDougall (1987)’s definition of (approximately) neutral surfaces (ANS). Now, if we use $B=\lambda$, in the inviscid case ${\bf F}=0$ and invoke Bernoulli theorem, $DB/Dt=0$, (50) predicts that

[TABLE]

Comparing this with (49) shows that

[TABLE]

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

3.4 Uniqueness of the inactive wind solution

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

[TABLE]

Since these two solutions both satisfy the continuity equation, a question that naturally arises is whether these may actually define the same velocity field despite being constructed from different fields. To examine this, let us recall that by definition ${\bf P}_{a}={\bf P}_{h}-\nabla B_{r}$, while $B_{a}=B_{h}-B_{r}$. It follows that

[TABLE]

which suggests that ${\bf v}_{ia}^{h}$ and ${\bf 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 $\boldsymbol{\omega}_{a}$. In other words, there must exist some scalar $\delta\lambda$ such that

[TABLE]

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

[TABLE]

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 indeed generally the case. It will be of interest to ascertain whether the vertical velocity component of ${\bf v}_{ia}^{a}$ is a better predictor of the actual vertical velocity field than that of ${\bf 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}=\boldsymbol{\omega}_{a}\cdot\nabla\gamma^{T}/\rho$ the PV variable constructed from it. As discussed in Saenz et al. (2015); Tailleux (2016, 2021); Tailleux and Wolf (2023), APE theory naturally comes with a generalised form of 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 Jackett and McDougall (1997) empirical neutral density variable $\gamma^{n})$, while also being mathematically and physically well defined. This motivated Tailleux (2016, 2021) to define thermodynamic neutral density $\gamma^{T}$ as

[TABLE]

where $f_{n}(p_{r})$ is a polynomial pressure correction empirically fitted to make $\gamma^{T}$ looks 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 generalisation of the concepts of density and PV to compressible seawater, then one may proceed similarly as in Needler (1985) 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

[TABLE]

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

[TABLE]

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, as it is one that is a priori as equally valid for an ideal fluid, for which it was originally developed by Needler (1985), as well as for a fully compressible ocean with a realistic nonlinear equation of state, as recently initiated by Ochoa et al. (2020). 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. Such a solution was previously derived by Gassmann (2014) in the context of a dry atmosphere and 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 Muller (1995) (that is itself related to the generalised Bernoulli theorem of Schär (1993)). 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 ${\bf 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 of $S$ and $\theta$ (for compressible seawater). We find that only if the ‘available’ form of Bernoulli function is used it is 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 state of minimum potential energy entering Lorenz APE theory Lorenz (1955). Indeed, only in this case is the vector ${\bf P}$ parallel to an approximately neutral density variable, namely the Lorenz reference density (LRD) discussed at length in Tailleux (2016, 2021); Tailleux and Wolf (2023).

Our third main result is that inactive wind solutions defined for different forms of Bernoulli function and vector ${\bf 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 small Rossby number, this is equivalent to say that different inactive wind solution differ primarily in their vertical velocity component. More generally, this result means that while it is in principle 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 against 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 towards a more rigorous and general theory of the oceanic 3D velocity field valid for a realist compressible ocean, as we hope to further demonstrate through concrete applications in future studies.

\funding

This research has been supported by the NERC-funded OUTCROP project (grant no. NE/R010536/1).

\dataavailability

The WOCE Global Ocean Climatology 1990-1998 (file ’wghc_params.nc’) used in this study is available at doi:10.25592/uhhfdm.8987. Enthalpy and other functions of state were estimated using the TEOS-10 library available at www.teos-10.org. Software to compute the analytic form of thermodynamic neutral density is available at doi:10.5281/zenodo.4957697. See also https://github.com/GammaTN for regular software updates and illustrative code examples.

Acknowledgements.

The author gratefully acknowledges comments from A. Colin de Verdière and Geoff Stanley on a original version of the paper, as well as constructive and supportive comments from three anonymous referees. \conflictsofinterestThe author declares no conflict of interest \abbreviationsAbbreviations The following abbreviations are used in this manuscript:

[TABLE]

\appendixtitlesno \appendixstart \reftitleReferences