## Abstract

The materiality and neutrality of neutral density and several forms of orthobaric density are calculated and compared using a simple idealization of the warm-sphere water mass properties of the Atlantic Ocean. Materiality is the value of the material derivative, expressed as a quasi-vertical velocity, following the motion of each of the variables: zero materiality denotes perfect conservation. Neutrality is the difference between the dip in the isopleth surfaces of the respective variables and the dip in the neutral planes. The materiality and neutrality of the neutral density of a water sample are composed of contributions from the following: (I) how closely the sample’s temperature and salinity lie in relation to the local reference *θ*–*S* relation, (II) the spatial variation of the reference *θ*–*S* relation, (III) the neutrality of the underlying reference neutral density surfaces, and (IV) irreversible exchanges of heat and salinity. Type II contributions dominate but have been neglected in previous assessments of neutral density properties. The materiality and neutrality of surfaces of simple orthobaric density, defined using a spatially uniform *θ*–*S* relation, have contributions analogous to types I and IV, but lack any of types II or III. Extending the concept of orthobaric density to permit spatial variation of the *θ*–*S* relation diminishes the type I contributions, but the effect is counterbalanced by the emergence of type II contributions. Discrete analogs of extended orthobaric density, based on regionally averaged *θ*–*S* relations matched at interregional boundaries, reveal a close analogy between the extended orthobaric density and the practical neutral density. Neutral density is not superior, even to simple orthobaric density, in terms of materiality or neutrality.

## 1. Introduction

It has been a century since the demonstration of the relation of ocean currents to density gradients (Helland-Hansen and Nansen 1909). The data from the *Meteor* expedition (Wüst 1935) showed the close alignment of core layers in the ocean, on which water properties like salinity, dissolved oxygen, etc., assumed extremal values, with strata of density corrected, however crudely, for adiabatic compression. [Usually, *σ _{t}* was used, which is merely the in situ density with its pressure dependence neglected; occasionally, the specific volume anomaly has been employed (Reid 1965).] Montgomery (1938) gave a cogent, qualitative explanation of why water properties should spread out along (corrected) density surfaces.

For a long time, potential density was accepted as the correct standard for comparison with water property core layers, but Lynn and Reid (1968) showed a sensitive and anomalous dependence of the topography of potential density surfaces on the choice of their reference pressure. The distortion of these surfaces complicates the relationship between potential density gradients and dynamically important pressure gradients (de Szoeke 2000). This is most dramatically illustrated by the well-known reversal of the vertical gradient of the potential density (referenced to zero pressure) in a stably stratified region of the deep Atlantic Ocean. The cause can be traced to the fact that the effective thermal expansion coefficient of potential density is that of the in situ density at the potential density’s reference pressure, thereby incorrectly representing thermal effects on buoyancy at pressures far from the reference pressure (McDougall 1987b). To mitigate this shortcoming, potential density has been computed over narrow pressure ranges (e.g., 0–500, 500–1500, 1500–2500 dbar, etc.), referenced to pressures within the respective ranges (*σ*_{0}, *σ*_{1}, *σ*_{2}, etc.) and matched at the transition pressures to form composite surfaces (Reid and Lynn 1971). A further elaboration allows for geographical variation by matching different values of potential density across the transition pressures in different regions. This procedure [sometimes called the patched potential density method (de Szoeke and Springer 2005)] has been used successfully (e.g., Reid 1994) to create global-scale maps of potential temperature and salinity, radio tracers such as tritium corrected for its radio decay, anthropogenic tracers like chlorofluorocarbons, or nutrients like nitrate and phosphate corrected for oxygen concentration [so-called NO and PO (Broecker 1981)].

The formulation of the patched potential density is not entirely satisfactory (de Szoeke and Springer 2005). A significant limitation is that it is defined discretely in both the vertical and horizontal directions. The vertical discretization was refined by de Szoeke and Springer (2005), who showed that decreasing the vertical pressure interval of the patches from the typical 1000 to 200 dbar, say, may render the matching procedure ill-defined because there is no longer a one-to-one correspondence between potential density surfaces across regional boundaries. Two major definitions of continuous corrected density variables have been proposed: neutral density (Jackett and McDougall 1997) and simple orthobaric density (de Szoeke et al. 2000). These two variables seem, at least superficially, to be very different. Part of the goal in this paper is to derive them within a common framework so that their relationship can be understood more easily. A second goal is to compare their performance in idealized settings that are analytically tractable. We propose three quantitative standards for comparison: 1) How material are the surfaces of the variable? 2) How closely do the gradients of the variable correspond to dynamically important density gradients? 3) How neutral are surfaces of the variable? The first two standards are motivated directly by early discoveries about the relationship between thermodynamic variables and motion, and the third is motivated by considerations of the preferred directions of mixing. These concepts will now be defined precisely.

Strict materiality is the quality of being material, or conservative:^{1} a fluid variable is material (conservative) if it does not change following the motion, or *Dχ*/*Dt* = 0. Quasi-materiality is the property of a variable changing only in response to irreversible, diabatic effects; that is, *Dχ*/*Dt* = *q̇ ^{χ}* ≡

**∇**· (𝗞

**∇**

*χ*), where 𝗞 is a diffusivity tensor for the variable. In this paper, by materiality we shall mean (at very slight risk of confusion) quasi-materiality. This is not to imply that irreversible contributions to the balance of

*χ*are negligible, only that reversible contributions to strict materiality are under consideration. Potential temperature and potential density (with a single reference pressure) are material variables. However, in situ density

*ρ*governed by the balance

is not material because of reversible adiabatic compression, apart from irreversible effects *q̇ ^{ρ}*.

Though not material, density (or specific volume *V* = *ρ*^{−1}) has the following crucial substitution property for the combination that occurs in the equations of motion of pressure gradient force and gravitational acceleration:

where *M*^{(V)} = *Vp* + *gz* is the Montgomery function (closely related to dynamic height) of density surfaces and **k**′^{(V)} = (**k** − **∇*** _{V}z*)/∂

*.*

_{V}z^{2}The exact replacement, as it were, of −

*ρ*

^{−1}

**∇**

*p*by the three-dimensional gradient of

*M*

^{(V)}, expressed relative to nonorthogonal slanting nonunit basis vectors

**i**,

**j**, and

**k**′

^{(V)}, and of −

*g*

**k**by

*p*

**k**′

^{(V)}in the vertical component, is of crucial consequence to the balance of forces regarded from isopycnal surfaces (de Szoeke 2000). The effective buoyant force acts in the direction of

**k**′

^{(V)}=

**∇**

*V*, that is, normal to isopycnals. When the left side of (1.2) is set equal to the Coriolis force, this substitution becomes the basis of the thermal wind calculation of the dynamic height and geostrophic currents:

*M*

^{(V)}= ∫

*p dV*(Sverdrup et al. 1942). Taking the vertical component of the curl of the momentum balance relative to isopycnal coordinates furnishes a potential vorticity based on the spacing of in situ isopycnals. The balance of this potential vorticity is marred, however, by very strong vortex stretching due to motions across in situ density surfaces arising from adiabatic compression.

The density of seawater is a function of potential temperature, salinity, and pressure. If there were a temperature–salinity relation unvarying in space or time in the oceans, density could be regarded as a function only of pressure and potential temperature *θ*: *ρ*^{−1} = *V*(*p*, *θ*). In that case the pressure gradient plus gravitational acceleration combination occurring in the momentum balance may be replaced by the following:

where *M*^{(θ)} = *gz* + *H*^{(θ)} is the Montgomery function, the sum of the gravitational potential and heat content, the latter defined by ; Π^{(θ)} = [∂*H*^{(θ)}/∂*θ*]* _{p}* is the Exner function; and

**k**′

^{(θ)}= (

**k**−

**∇**

*)/∂*

_{θ}z*=*

_{θ}z**∇**

*θ*. This has obvious parallels to (1.2), with similar interpretations. The potential vorticity based on

*θ*-surface (isentrope) spacing ∂

*is far more attractive than the*

_{θ}z*ρ*-surface potential vorticity because the potential temperature is material and the dientropic vortex stretching is much reduced or negligible (Eliassen and Kleinschmidt 1954). Also, the potential temperature gradient

**∇**

*θ*is parallel to the dianeutral vector,

(about which more will be developed below), when there is an unvarying temperature–salinity relation. Isentropes in this situation are neutral surfaces. This property is closely linked to the material property of potential temperature.

The nexus among the three properties of potential temperature—dynamical substitution (1.3), materiality, and neutrality—is broken when the temperature–salinity relation is not unvarying but dependent on geographic location. In such a case the dianeutral vector **d** cannot be parallel everywhere to any scalar gradient, let alone **∇***θ*, though it may be taken to define at any point a neutral plane to which it is perpendicular. An early hope was that the totality of neutral planes would define a generalization of neutral surfaces, but this does not happen (McDougall 1987a; Davis 1994). The formation of well-defined surfaces from neutral planes requires that their helicity vanish everywhere; that is, **d** · **∇** × **d** = 0. This amounts to demanding an unvarying temperature–salinity relation, which brings us back to isentropes, closing the logical circle! Failing this, global measures of the closeness to neutral planes of the isopleth surfaces of a given variable can be devised—which could be called the surfaces’ neutrality. For example, Eden and Willebrand (1999) considered several global indices of the mean-square difference, variously weighted, between slopes of candidate surfaces and the neutral plane slope. A variable that labels such candidate surfaces may be called a neutral density of the first kind.

Simple orthobaric density *σ*^{⊕}—a function only of density and pressure, and called simple to distinguish it from extended orthobaric density, considered below (section 5)—is calculated by assuming a virtual compressibility *ρ*^{−1}*c*_{0}(*ρ*, *p*)^{−2}—also a function only of pressure and density—to correct the in situ density for compression: *ϕ dσ*^{⊕} = *dρ* − *c*_{0}(*ρ*, *p*)^{−2}*dp*, where *ϕ*(*ρ*, *p*) is an integrating factor. De Szoeke et al. (2000) fitted a function of *p* and *ρ* to the global field of sound speed in the ocean to furnish *c*_{0}(*ρ*, *p*). This was shown to be equivalent to fitting a global average *θ*–*S* relation. Because it depends only on the in situ density and pressure, the simple orthobaric density has the dynamical property (1.3), with *σ*^{⊕} substituted for *θ*. The material derivative of *σ*^{⊕} may be calculated (de Szoeke and Springer 2005). The simple orthobaric density is not material unless the true sound speed *c* is identical with the hypothetical *c*_{0} everywhere; that is, unless the assumed *θ*–*S* relation actually applies everywhere. The difference *c*_{0} − *c* is a measure of the materiality of the orthobaric density. A normalized index *ψ* of materiality, called the buoyancy gain factor, was devised for simple orthobaric density. This factor is so called because it gives the ratio of the apparent buoyancy frequency squared, based on the orthobaric density, to the true buoyancy frequency squared. Its difference from 1 is proportional to *c*_{0} − *c* and is, thereby, a measure of the ocean’s deviation from the hypothetical universal *θ*–*S* relation. One of the tasks of this paper is to quantify the neutrality of the simple orthobaric density (section 4). Another task is to explore how using spatially variable *θ*–*S* relations to define the extended orthobaric density affects materiality and neutrality (section 5).

Jackett and McDougall (1997) described how to compute a variable *γ*, which they called neutral density, constructed in such a way that its isopleth surfaces^{3} are closely tangent to neutral planes everywhere. Given, say, a synoptic hydrographic section, the method calculates apparent displacements of water samples on the section from a standard global reference dataset, such as the Levitus climatology, and associates neutral density values preassigned to the reference data to the synoptic samples. One of the goals of this paper is to clarify the definition of this neutral density, designated here as the second kind, and, in particular, to distinguish it from neutral density of the first kind. McDougall and Jackett (2005a) proposed measures of materiality and neutrality for neutral density (section 2), modeled on de Szoeke et al.’s (2000) work with orthobaric density. However, we will show that in their calculation of neutral-density materiality a consequential term was unjustifiably neglected, and their conclusions about the relative materiality of neutral density and orthobaric density (McDougall and Jackett 2005b) are suspect. In addition, we quantify the neutrality of the neutral density for the first time and derive an expression for the additional forces that arise on neutral density surfaces, which do not possess the substitution property (1.3) (see the appendix).

The plan of this paper is as follows. In section 2 we briefly review how the neutral density of a water sample is calculated relative to a reference global hydrographic dataset. In the appendix, we calculate the material derivative of neutral density. We also derive a formula for the difference between the slopes and orientations of neutral density surfaces and neutral planes (neutrality), and a generalization of the dynamical substitution property (1.3) for neutral density. In section 3 we take as the reference dataset for neutral density a simple model of the Atlantic that represents the ocean’s meridional variation of the *θ*–*S* relation, though constructed in such a way that in situ isopycnals are level and the reference ocean is motionless. Supposing the “real” ocean to have been vertically displaced from the reference ocean by meridionally varying amounts, we calculate, by methods analogous to Jackett and McDougall’s (1997), the distribution of neutral density and its material derivative. In section 4 we calculate the simple orthobaric density, based on a mean *θ*–*S* relation of the reference dataset, of the displaced real ocean, as well as its material derivative and neutrality, which we compare to those of neutral density. We also calculate the dual-hemisphere orthobaric density (de Szoeke and Springer 2005), where different mean northern and southern *θ*–*S* relations representative of the respective hemispheres are employed. We show how the two different simple orthobaric densities in the hemispheres may be patched together. In an extension of these ideas, by subdividing the ocean into more and more provinces, each with its distinctive mean *θ*–*S* relation, and eventually passing to the obvious continuous limit, we obtain what we have called the extended orthobaric density (section 5), whose properties can be compared to those of neutral density. A discussion and summary conclude the paper (section 6).

## 2. Neutral density

### a. Preliminaries

In what follows it is useful to have the following definitions and concepts at hand. The equation of state *ρ*^{−1} = *V*(*p*, *θ*, *S*) is a function of pressure, potential temperature, and salinity (Jackett et al. 2006). The dianeutral vector defined by (1.4) may also be written

(McDougall 1987a). Its vertical component is **k** · **d** = *N* ^{2}/*g* (where **k** is a unit vertical vector, *g* is the earth’s gravitational acceleration, and *N* is the buoyancy frequency) and *a*, *b* are the thermal expansion and haline contraction coefficients. The dip of a neutral plane, or neutral dip, is the two-dimensional vector whose direction in the horizontal gives the orientation of the plane, and whose magnitude is the tangent of the angle between the neutral plane and the horizontal, namely,

[This usage is common in describing bedding planes in physical geology (Holmes 1945).] The dianeutral vector is integrable, that is, **d** = *ϕ***∇***ν*, so that **h** = −**∇*** _{z}ν*/∂

*=*

_{z}ν**∇**

*, if and only if the helicity vanishes:*

_{ν}zwhere *μ*_{*} = ∂ln(*a*/*b*)/∂*p*. The latter has typical reciprocal values for seawater of (McDougall 1987b; Akitomo 1999; Eden and Willebrand 1999). Condition (2.3) is satisfied if *S* = *S*(*θ*, *p*). In that case, surfaces of constant *ν* are everywhere tangential to neutral planes (perpendicular to **d**) and may be properly called neutral surfaces. The dynamical substitution (1.3) is valid, with *ν* substituted for *θ*. (Under the slightly more restrictive condition *S* = *S*(*θ*), *ν* may be taken to be identical to *θ*, as was assumed in the introduction.)

Various practical kinds of neutral density have been devised even when (2.3) is not satisfied. But the tangent planes of their surfaces cannot perforce be perpendicular everywhere to the dianeutral vector, so estimates of their closeness to this “ideal” ought to be formed. This is one of the tasks undertaken in the following.

### b. Neutral density of the first kind

Suppose we call **∇*** _{γ}z* the dip of a surface

*z*=

*z*(

*γ*,

**x**), as labeled by the constant parameter

*γ*. The vector difference between this dip and the neutral dip,

is called the neutrality of the surface. Its magnitude |**∇*** _{γ}z* −

**h**| gives the angle between the surface and the neutral plane. When the neutrality is zero, surfaces of constant

*γ*are said to be neutral. Even if not, we shall call these surfaces neutral density surfaces of the first kind and call

*γ*the neutral density,

^{4}when

**∇**

*≈*

_{γ}z**h**, to some acceptable degree of approximation. Several ways of approximating this equality have been described. One way is to integrate

where *s* is the distance along a horizontal path **x**(*s*), to form a three-dimensional neutral trajectory **x**(*s*), *z*(*γ*, *s*). By repeating this over a network of such paths chosen to cover the entire ocean, and “connecting” paths sharing the same value of *γ*, McDougall (1987a) suggested that “approximate surfaces” may be formed. There are unavoidable mathematical difficulties with this approach when the helicity of the dianeutral vector does not vanish (Davis 1994). However, a realizable approach is to choose *z*(*γ*, **x**) to minimize the neutrality globally, that is, to minimize

where *V* is the total volume of the ocean and **W** is a suitable weight function. Eden and Willebrand (1999) considered several similar minimization principles. Neutral density surfaces of the first kind exist even when the neutrality **∇*** _{γ}z* −

**h**differs from zero. They must be regarded, however, as one step removed from neutral planes.

### c. Neutral density of the second kind

It is impractical to recalculate the global neutral density field whenever a new hydrographic section becomes available. Accordingly, Jackett and McDougall (1997) described a two-stage procedure for assigning the value of a second kind of practical neutral density to a sample of water collected at a certain horizontal position and ambient pressure.

- (i) The first stage is to prepare beforehand a set of reference neutral density surfaces of the first kind from an archive of hydrographic data, such as the Levitus climatology, from which averaged vertical reference profiles of potential temperature and salinity,
*θ*_{ref}(*p*,**x**),*S*_{ref}(*p*,**x**),^{5}at any position**x**in the world’s oceans have been compiled. (One might call this the reference ocean.) The reference surfaces are specified by the function - (ii) Next, suppose a water sample has been collected with the values
*θ*,*S*at ambient pressure*p*, position**x**, and time*t*, from the present ocean (so called to distinguish it from the reference ocean). The target pressure level*p*is the pressure at which the reference ocean’s potential density matches the present sample’s potential density, each referenced to the weighted-average^{r}*p*=*λ̃p*+*λp*[with^{r}*λ̃*+*λ*= 1; Jackett and McDougall (1997) invariably use*λ*=*λ̃*= 0.5 for the weights]. This prescription may be written as and is called the neutral property (Jackett and McDougall 1997). It furnishes*p*as a function of^{r}*θ*,*S*,*p*, and position**x**. The last step is to obtain the neutral density of the second kind of the water sample by substituting*p*into the predetermined function^{r}*γ*_{ref}(*p*,^{r}**x**).

Graphically, this identifies the reference neutral density surface that goes through the target level *p ^{r}* at position

**x**in the reference ocean, and associates its neutral density (first kind) with the present water sample. This second kind of neutral density is two steps removed from neutral planes: first by the association the neutral property [Eq. (2.8)] makes between the present ocean and the reference ocean, and second by the relation of the reference neutral density to neutral planes.

### d. Material rate of change: Cross-neutral mass flow

An important property of neutral density of the second kind can be obtained by taking the material derivatives^{6} of both sides of Eqs. (2.8) and (2.7) and eliminating *Dp ^{r}*/

*Dt*to furnish an equation for its material rate of change,

*γ̇*≡

*Dγ*/

*Dt*, in the present ocean. (Details are given in the appendix.) We call the value of

*γ̇*at a point in space and time the materiality; if the materiality were zero everywhere,

*γ*would be perfectly material, or conservative. Materiality gives the rate of flow

*e*(m s

^{γ}^{−1}, positive when upward) across neutral density surfaces. This is obtained in the appendix by Eq. (A.2), which is reproduced here:

In this equation, *ψ̂ ^{γ}* = 1 −

*μ*[

*ϑ*−

*ϑ*

_{ref}(

*p*,

^{r}**x**)]

*D*

^{−1}is the neutral-density buoyancy gain factor;

^{7 }

*ϑ*,

*s*are normalized forms of

*θ*,

*S*(slightly nonlinearly transformed in the case of

*ϑ*; see the appendix);

*μ*=

*V*/

_{θp}*V*is the thermobaric parameter (McDougall 1987b);

_{θ}*D*is a factor given by (A4) in the appendix; and Δ

*σ*∂

*p*/∂

*γ*is the vertical thickness (in pressure units, i.e., decibars) of an increment Δ

*σ*of neutral density. The right side of (2.9) uses a simplified equation of state given by Eqs. (A.1a) and (A.1b) in the appendix.

One may distinguish four types of contributions to the material rate of change, or the diapycnal material flow, arising from the following sources:

(I) the sum of the apparent vertical motion in the present ocean,

*w*= (∂_{a}+_{t}**u**·**∇**), as though the constant-_{γ}p*γ*surface were material, and in the reference ocean,*w*=_{a}^{r}**u**·**∇**(Fig. 1), each weighted by_{γ}p^{r}*λ̃*or*λ*, respectively, and attenuated by*ψ̂*− 1, which is a measure of the departure of the water sample’s^{γ}*ϑ*,*s*from the local reference*θ*–*S*relation;(II) differences between ambient pressure in the present data and target pressure in the reference dataset,

*p*−*p*, combined with advection of the reference temperature on neutral density surfaces (if the reference temperature changes on the reference neutral density surfaces, it is because of the spatial variation of the reference ocean’s^{r}*θ*–*S*relation);(III) reference neutrality, that is, differences between the neutral density surface dip

**∇**_{γ}z_{ref}= −**∇**/(_{γ}p^{r}*ρ*_{ref}*g*) and the true neutral dip**h**_{ref}(see the appendix) in the reference dataset; and(IV) the irreversible rates of change of temperature and salinity that constitute the buoyancy rate of change

*q̇*^{(p)}[see Eq. (A5) in the appendix].

### e. Neutrality

Another result established in the appendix is an equation for the neutrality, the difference in the present ocean between the dip **∇*** _{γ}z* of the neutral density (of the second kind) and the neutral dip

**h**. This difference is given by Eq. (A.9):

One may identify contributions analogous to types I, II, and III for materiality. Note the distinction between the neutralities of the present ocean (left-hand side) and the reference ocean (term III).

McDougall and Jackett (2005b) derived an expression similar to (2.9). Jackett and McDougall (1997) argued that, given their method of preparing the reference dataset, the type III term in (2.9) [and by inference in (2.10)] is small. McDougall and Jackett (2005a) conducted a global census of *ψ̂ ^{γ}* from the world hydrographic dataset showing that it is very close to 1 (

*ϑ*−

*ϑ*

_{ref}≈ 0) and concluded that the type I contribution to the materiality (and, one may infer, neutrality) of neutral density due to deviations from the local

*θ*–

*S*relation is negligibly small. Further, following from a claim that

*μ*(

*ϑ*−

*ϑ*

_{ref})

**u**·

**∇**

*(type I) and*

_{γ}p^{r}*μ*(

*p*−

*p*)

^{r}**u**·

**∇**

_{γ}ϑ_{ref}(type II) were typically of the same order, they asserted that diapycnal flow contributions from spatial gradients of the reference

*θ*–

*S*relation were likewise negligible. This claim is very dubious. In section 3 we shall present an illustration clarifying that type II contributions to neutrality and materiality cannot be neglected.

### f. Neutral density in the displaced reference ocean

The treatment in the appendix is quite general with regard to the reference dataset and the present ocean. A special case is of particular interest. Suppose the present ocean is obtained from the reference ocean by displacing the latter vertically by an increment Δ*p*, possibly itself a function of *p*, **x**, and *t*:

When this is substituted into the right side of (2.8), an obvious solution for target pressure is

## 3. Properties of neutral density in an idealized ocean

To evaluate the performance of neutral density (and, in section 4, orthobaric density) with regard to materiality and neutrality, it is useful to consider a simplified model of the Atlantic Ocean’s meridional temperature and salinity structure that illustrates its salient hydrographic features. This model will be taken as the reference dataset for the neutral density algorithm. Its motivation is the spatial dependence of the Atlantic’s *θ*–*S* relation, which was identified in the preceding sections as the cause of nonzero helicity and the nonexistence of neutral surfaces. We will suppose the present ocean to be a vertically displaced form of the reference ocean from which we can calculate the neutral density distribution and its materiality and neutrality.

### a. A simple model of the warm-water Atlantic Ocean

The *θ*–*S* relation of the Atlantic on the horizontal 500-dbar pressure surface is shown in Fig. 2 in terms of the modified variables *ϑ*, *s* [Eq. (A1b) in the appendix]. The diagram shows the main southern and northern lines of the tropical to temperate upper ocean, with a constant-density bridge linking the main lines. The nonlinear transformation of potential temperature *θ* has straightened a slight curvature of the lines. On any pressure surface shallower than 1500 dbar [only 500 dbar is shown here, but see de Szoeke and Springer (2005) for more pressure surfaces], the main lines are practically identical, while the bridge lies at different densities. The two fitted lines in Fig. 2,

may be considered idealized representations of the main *θ*–*S* lines at temperatures above 4°C (*ϑ* ≥ −0.11 × 10^{−3}). More elaborate representations of the *θ*–*S* lines than (3.1), taking account of the convergence of the northern and southern lines at colder temperatures, may readily be considered, though in the interests of brevity we have foregone such a treatment here. It may be more conventional to show the spatial variation of *θ*–*S* relations as in Fig. 3, where each line of data represents a meridional location. This illustrates the evolution along a meridional section from the southern to the northern extremes idealized in Eq. (3.1). The near-constant density bridge structure at each pressure, exemplified in Fig. 2, unfolds to form the meridionally dependent *θ*–*S* profiles of Fig. 3. The fit of parallel, straight lines to the *ϑ*–*s* relation at each latitude for waters warmer than 4°C or shallower than 1500 dbar seems quite accurate. The lines are biased to the saltier side on passing from 20°S to 20°N, but are unvarying with latitude south and north, respectively, of these limits (overlooking some anomalously salty water from the Mediterranean at northern midlatitudes). A linear fit to this bias between 20°S and 20°N may be oversimplified but not unrepresentative.

Thus, Figs. 2 and 3 motivate the following idealizations of the temperature and salinity fields of the warmer layers of the subtropical to tropical Atlantic:

for |*y*/*L*| ≤ 1, where *y* = ±*L* are the northern and southern geographical boundaries of the transequatorial Atlantic bridge region. For *y*/*L* > +1 (*y*/*L* < −1), we replace *y*/*L* in (3.2a) by +1 (−1). We take 2*L* ≈ 3000 km and *n* = 1 − *R*^{−1} ≈ 0.5. (The origin, *y* = 0, need not be the equator.) The *ϑ*_{ref} and *s*_{ref} fields are chosen to give a horizontally uniform profile of the normalized specific volume anomaly,

here assumed linear with depth with constant coefficient *σ*′* _{B}*. We call this the quiescent reference ocean (Fig. 4). One might as well assign the

*δ*values of level pressure surfaces (or rather their negatives) as the neutral density labels of the reference ocean, that is, as in Eq. (2.7):

^{8}

If *δ*_{ref} is eliminated from (3.2a) and (3.2b), a meridionally variable reference *ϑ*–*s* relation emerges:

(cf. Fig. 3). The *ϑ*–*s* relations at successive positions from *y* = −*L* to *L* are straight lines of slope *R*, shifting to saltier water. The total shift is 2*s*_{1}. South of *y* = −*L* (or north of *y* = *L*) the *ϑ*–*s* relation is taken to be spatially uniform, as given by Eq. (3.1). The end-point northern and southern *ϑ*–*s* relations, though very simple, represent the real ocean very well, as does the mapping of the intermediate *ϑ*–*s* relations onto the constant-density bridges (Fig. 2). More schematic than strictly accurate is the assumption of a meridionally linear transition between the end points (Fig. 3).

A consequence of the density-compensated bridge structure [Eq. (3.2)] is that in situ isopycnals and neutral planes (McDougall 1987a) are horizontal in the reference dataset. The dianeutral vector [Eq. (2.1)] is everywhere vertical,

hence, the reference helicity is trivially zero: **d**_{ref} · **∇** × **d**_{ref} = 0. The reference neutrality is also zero; the reference state exactly (and trivially) satisfies −**∇*** _{γ}p^{r}*/(

*ρ*

_{ref}

*g*) =

**h**

_{ref}= 0, so that the type III terms that occur in (2.9) and (2.10) are absent. Reference-neutral surfaces are level.

From Eq. (3.2) describing the quiescent reference state, one obtains

where **j** is the meridional unit vector. [When approximations are given in (3.3), the contributions from the small parameter *μ* have been neglected.]

### b. Idealized neutral density in the displaced reference ocean

Let the place of the hydrographic archive used as the reference data for the neutral density algorithm (2.8) be taken by the two-dimensional quiescent reference state [Eq. (3.2)]. Suppose that the present ocean is the reference ocean displaced by Δ*p*(*p*, **x**), as in Eqs. (2.11). In general, then, helicity will not be zero and neutral surfaces will not exist. However, neutral density and its isopleth surfaces are well defined and can be obtained by inserting Eq. (2.11) on the left of (2.8). As remarked above, the target pressure *p ^{r}* is then given by (2.12), and the neutral density by

Neutral density surfaces are described by *p ^{r}* = const. Interesting properties of neutral density, such as materiality and neutrality, may be estimated. The buoyancy gain factor for the displaced ocean is

*ψ̂*= 1 (because

^{γ}*ϑ*−

*ϑ*

_{ref}(

*p*,

^{r}**x**) = 0), and type I terms in (2.9) and (2.10) vanish. In a pressure range centered around 500 dbar (for

*p*<

*p*

_{0}∼ 1000 dbar, say), suppose the displacement is linear in

*x*and

*y*:

The parameter ɛ in (3.4) may be so small as to be negligible in Δ*p* in (3.5), but it is essential in providing a nonzero meridional velocity *υ* = *f* ^{−1}*σ*′* _{B}*(

*p*−

*p*

_{0})Δ

_{0}ɛ/2

*L*to make (3.5) nonzero and nontrivial. Fundamentally, the zonal structure of (3.4) provides three-dimensional structure in the displaced ocean to furnish nonzero helicity. The meridional dependence of this is shown in Fig. 5a for

*λ̃*= 0.5, normalized by a parameter including the meridional velocity

*υ*.

The neutrality of neutral density is given by Eq. (2.10). Applying this formula to the displaced reference ocean, in which only the type II term on the right survives, one obtains, using the approximations of (3.3),

This difference, or mismatch, is shown in Fig. 5b. Curiously, even if ɛ ≡ 0, despite the displaced ocean’s being perfectly two-dimensional with zero helicity (so that neutral surfaces are well defined), the neutrality of neutral density surfaces does not vanish. An estimate of the dimensionless parameter appearing in both (3.5) and (3.6) is

This gives *e ^{γ}* ∼ 10

^{−5 }

*υ*∼ 10

^{−7}–10

^{−6}m s

^{−1}for

*υ*∼ 10

^{−2}–10

^{−1}m s

^{−1}. In comparison to typical Ekman pumping magnitudes near the ocean surface (∼10

^{−6}m s

^{−1}), such dianeutral material flows are not inconsiderable.

## 4. Properties of orthobaric density in an idealized ocean

The idealized representation in section 3 of the Atlantic Ocean illustrates well the central irremediable feature of the ocean’s density field that gives nonzero helicity: the spatial variation of the *θ*–*S* relation. We now turn to calculating orthobaric density functions predicated on *θ*–*S* relations from various locations in the reference ocean. These relations and density functions are summarized in Table 1. We also calculate their buoyancy gain factors (Table 2), cross-isopleth mass flows (materiality), and neutrality (difference of the surfaces’ **∇***z* from neutral dip **h**), and compare them with those for neutral density.

### a. Simple orthobaric density

An orthobaric density function for the entire idealized Atlantic may be defined in the following way. Take the reference *ϑ*–*s* relation (3.2e) at some central location, say *y* = 0, that is, *s*_{ref} = *R*^{−1}*ϑ*_{ref}. By (3.2a), *ϑ*_{ref}|_{y=0} is related to a normalized specific volume anomaly *δ* [Eq. (A1a) in the appendix] and pressure *p* by

By substituting this into the adiabatic compressibility anomaly (∂*δ*/∂*p*)_{ϑ,s} = *μϑ*, one obtains the virtual compressibility anomaly, which may be used to pose the following differential equation:

This furnishes the integral −*δ*/(1 + *μp*/*n*) = const. One may take the value of this “constant,” defining a surface in *δ*, *p* space, at a reference pressure, say *p* = 0, to be *σ*^{⊕}. In this way one obtains a simple “global” orthobaric density^{9} as a function of *δ*, *p*:

The material rate of change of *σ*^{⊕} may be calculated. It can be written as

(de Szoeke et al. 2000; de Szoeke and Springer 2005). Here, Δ*σ*∂*p*/∂*σ*^{⊕} is the vertical thickness (dbar) of an increment Δ*σ* of orthobaric density. In addition, *e*^{⊕} is the net flow (m s^{−1}) across a *σ*^{⊕} isopycnal and consists of two parts. The first is an adiabatic part—a residual of the compressibility effect—given by the apparent motion (∂_{t} + **u** · **∇**)_{σ⊕}*p* of orthobaric isopycnals, attenuated by *ψ*^{⊕} − 1, where

is the buoyancy gain factor, in which

Only at the location from which the reference *ϑ*–*s* relation has been taken (*y* = 0) will *ϑ* = *ϑ*_{ref} in the displaced ocean, and *ψ*^{⊕} = 1, reducing term I in (4.2) to zero. The remaining part of (4.2), labeled IV, is a term proportional to *q̇*^{(p)}, the sum of all diabatic, irreversible influences on the negative buoyancy [Eq. (A5), see the appendix). We have assigned the terms in (4.2) type designations based on their appearance analogous to Eq. (2.9). Note that terms of types II or III do not appear.

The buoyancy gain factor of simple orthobaric density for the displaced reference ocean is listed in Table 2 and shown in Fig. 6. Note the linear decrease with *y* of *ψ*^{⊕} − 1. A similar decrease is evident in an Atlantic meridional section of *ψ*^{⊕} − 1 (de Szoeke et al. 2000, their Fig. 11). (Since the reference *θ*–*S* relation in the latter calculation was taken from southern waters, there is a positive bias in the real-world section compared to the idealized section.) Unlike the type I term in the neutral density materiality, Eq. (2.9), the contribution of the term of which this is a factor to material flow across orthobaric isopycnals does not vanish. Consider the steady-state estimate of the apparent vertical motion, (∂_{t} + **u** · **∇**)_{σ⊕}*p* ≈ *υ*(∂*p*/∂*y*)_{σ⊕} ≈ *υ*Δ_{0}/2*L*. The contribution from the type I term of (4.2) to diapycnal flow is

The meridional dependence of this, normalized by *μs*_{1}Δ_{0}/(*ρgnσ*′* _{B}L*)

*υ*, is shown in Fig. 5a; also shown is the meridional average value of |

*e*

_{I}

^{⊕}| (similarly normalized). One may compare

*e*

_{I}

^{⊕}to , given by (3.5). Though their meridional dependences are offset one from the other, their magnitudes and parameter dependences are similar. The latter, , was neglected by McDougall and Jackett (2005a). From (2.9) and (4.2),

*e*

_{I}

^{γ}∝

*ψ*− 1 and

^{γ}*e*

_{I}

^{⊕}∝

*ψ*

^{⊕}− 1. McDougall and Jackett (2005a) show persuasively that

*ψ*− 1 ≪

^{γ}*ψ*

^{⊕}− 1, so that indeed

*e*

_{I}

^{γ}≪

*e*

_{I}

^{⊕}. In our illustration, in fact,

*e*

_{I}

^{γ}≡ 0,

*ψ*− 1 ≡ 0. But in compensation it must be recognized that is not negligible.

^{γ}The dip of an orthobaric density surface −**∇**_{σ⊕}*p*/(*ρg*) differs from the neutral plane dip **h** [Eq. (2.2)] by

(de Szoeke et al. 2000). The last approximate equality follows from the first row entry of Table 2, and has a magnitude similar to the neutrality of neutral density surfaces [Eq. (3.6)]. See Fig. 5b.

The panels in Fig. 5 illustrate a major result of this paper. They show typical estimates of diapycnal flow (materiality) and the neutrality of neutral density and orthobaric density. The dianeutral material flow *e ^{γ}* must be expected to have magnitudes similar to the orthobaric diapycnal flow

*e*

^{⊕}. Likewise, the neutralities of neutral density and orthobaric density are similar.

### b. Regional orthobaric density functions

Orthobaric density functions can be defined specific to various regions, each predicated on a reference *ϑ*–*s* relation typical of its region. To construct global corrected density surfaces, we must specify how to match the various density functions across their interregional boundaries (de Szoeke and Springer 2005). We will illustrate this difficulty with some examples.

#### 1) Northern orthobaric density

An orthobaric density might be predicated on the *ϑ*–*s* relation (3.2e) at *y* = *L*, *s*_{ref} = *R*^{−1 }*ϑ*_{ref} + *s*_{1}. This choice is important because it characterizes the water mass properties of much of the northern North Atlantic, and the northern portion of the tropical bridge region. If such a choice is pursued along the preceding lines, one arrives at the specification of what might be called the northern orthobaric density function *σ ^{N}*, given in the second row of Table 1. [Any function of (−

*δ*−

*s*

_{1})/(1 +

*μp*/

*n*) might be chosen for

*σ*. We have chosen a particular linear function, for a reason that will become clear presently.] The material derivative of

^{N}*σ*gives an equation exactly like (4.2), involving a buoyancy gain factor

^{N}*ψ*, defined like (4.3), with ⊕ superscripts replaced by

^{N}*N*superscripts, and with

*ϑ*

_{ref}|

_{y=L}(

*σ*) given in Table 1. The buoyancy gain factor

^{N}*ψ*and the integrating factor

^{N}*ϕ̃*for this northern orthobaric density are listed in Table 2 for the displaced reference ocean;

^{N}*ψ*is shown in Fig. 6. Buoyancy gain factors on a meridional Atlantic section based on a northern (and a southern)

^{N}*θ*–

*S*relation are shown in de Szoeke and Springer’s (2005) Fig. 13.

#### 2) Southern orthobaric density

All this might be repeated to develop another orthobaric density function *σ ^{S}* based on the

*ϑ*–

*s*relation in (3.2e) at the southern end of the tropical bridge,

*y*= −

*L*; that is,

*s*

_{ref}=

*R*

^{−1}

*ϑ*

_{ref}−

*s*

_{1}. This will yield the entry in the third row of Table 1, a linear function of (−

*δ*+

*s*

_{1})/(1 +

*μp*/

*n*). Again, the material derivative of

*σ*, analogous to (4.2), may be obtained. [See also de Szoeke and Springer’s (2005) Fig. 13.]

^{S}#### 3) Dual-hemisphere orthobaric density

The constant *m* factors specifying *σ ^{N}*,

*σ*in Table 1 have been chosen so that when the undisplaced reference bridge specific volume anomaly (SVA) profile applies,

^{S}*δ*= −

*σ*′

*, the two are numerically identical. Thus, the overall dual-hemisphere quasi-orthobaric density function defined by*

_{B}pwill be continuous across *y* = 0 if the said density profile pertains there. However, if the density profile differs from this standard, say, *δ* ≈ −*σ*′* _{B}*(

*p*+ Δ

*p*), as in the displaced ocean of section 2, the difference between the two expressions of (4.6) on either side of

*y*= 0 is approximately (using entries from Table 1)

(In the approximate equality, higher-order terms in *μ* are neglected.) This gives a discontinuous displacement for a *σ ^{O}* isopycnal across

*y*= 0 of

(deeper on the north side, Δ*z* < 0 when displacement from the standard profile is upward, or Δ*p* > 0; see Fig. 6). This discontinuity is a potential site of material leakage. Given meridional velocity *υ*, there is volume flux through this discontinuity of 2*s*_{1}*μ*/(*ρgnσ*′* _{B}*)(

*v*Δ

*p*)

_{y=0}(m

^{2}s

^{−1}, positive when toward lighter water). Divided by 2

*L*, this gives an equivalent cross-isopycnal mass flow, averaged over the bridge region, of

## 5. Extended orthobaric density

The bridge region need not be split into only two parts, northern and southern, as for the dual-hemisphere density function. Suppose the bridge region is subdivided into a number *N* of intervals |*y* − *y*_{j}| < *L*/(*N* − 1), where *y*_{j} = *L*(2*j* − *N* − 1)/(*N* − 1) for *j* = 1, … *N*. The treatment outlined above for the dual-hemisphere density function may be applied, with obvious modifications, to each such subinterval, characterized by the *θ*–*S* relation given by (3.2e) at its center, *y* = *y*_{j}. This leads to the specification in each subinterval of an orthobaric density function *σ*^{(yj)} listed in the fourth row of Table 1. Each such density function has an associated buoyancy gain factor *ψ _{j}* (Table 2). The boundaries between the subintervals correspond to

*y*= 0 in the dual-hemisphere orthobaric density. The

*ψ*− 1 factors determining the quasi-materiality of a local orthobaric density, defined by reference to the local

_{j}*θ*–

*S*relation, are reduced by a factor of 1/(

*N*−1). The discontinuities, analogous to (4.7), then occurring at the edges of subintervals,

*y*

_{j+1/2}=

*L*(2

*j*−

*N*)/(

*N*− 1), are similarly reduced,

but are *N*−1 times more numerous (Fig. 7). Given a meridional flow *υ*, there will be a mass transport through each such gap, which, averaged over the spacing between gaps, 2*L*/(*N*−1), may be thought of as an equivalent diapycnal flow,

distributed over the interval |*y* − *y*_{j+1/2}| < *L*/(*N* − 1). These flows are shown in Fig. 8 for *N* = 2 (dual-hemisphere orthobaric density), 4, and 10; also shown are the (meridionally unvarying) averages of |*e*_{I}^{O}| ∼ *O*[1/(*N* − 1)] in each subinterval. In the limit *N* → ∞, tends to the continuous line shown in Fig. 8.

An alternative view is gained by starting from the continuous limit of the arbitrarily centered extended orthobaric density function (fourth row, Table 1):

The reference temperature field (3.2a) may be expressed in terms of this entity as

Otherwise, in the state displaced by Δ*p* from the quiescent,

This approximation will be used below.

### a. Material derivative

The material rate of change of (5.2) may be written

where

Equation (5.4) gives the rate at which matter crosses a *σ ^{O}* isopycnal. The terms labeled I and IV, a residual compressibility effect and a diabatic contribution, respectively, are analogous to similar terms in (4.2). For the displaced reference ocean, it can be easily shown that

*ϑ*−

*ϑ*

_{ref}(

*σ*;

^{O}*y*) = 0, so that

*ψ*= 1, and the first term in (5.4) vanishes. The contribution of term II may be written

^{O}It is proportional to the meridional motion *υ*, the meridional temperature gradient on *σ ^{O}* surfaces ∂

*ϑ*

_{ref}(

*σ*;

^{O}*y*)/∂

*y*≈

*s*

_{1}/

*nL*, the thermobaric coefficient

*μ*, the reciprocal buoyancy frequency squared ∂

*p*/∂

*σ*≈

^{O}ϕ̃^{O}ψ^{O}*ρg*

^{2}/

*N*

^{2}, and the displacement Δ

*p*from the quiescent state. It coincides with the continuous limit of (5.1), as shown in Fig. 8. The neutral density diapycnal flow [Eq. (3.5); see Fig. 8c] differs only in the factor

*λ̃*.

Reviewing the progression from simple global orthobaric density to dual-hemisphere orthobaric density, to *N*-fold discrete extended orthobaric density (EOB), to continuous EOB, one sees that type I diapycnal flow, the sole reversible contributor to diapycnal flow across simple orthobaric density surfaces, is augmented by a type II contribution when dual-hemisphere EOB is considered, of about the same magnitude. As the number *N* of subintervals is increased, the type II contributions supplant the type I contributions, which dwindle inversely with *N*. In the continuous limit, only type II contributions to the EOB diapycnal flow remain, resembling strongly, other than the factor *λ̃*, the diapycnal flow for neutral density. McDougall and Jackett (2005a), while confirming the smallness of the type I contributions, completely neglected the type II contributions to neutral density materiality. Their conclusion that neutral density is superior to orthobaric density in its material properties rests on this fallacy.

### b. Neutrality

The difference between the dips of extended orthobaric density surfaces and neutral planes can be calculated. Take the quasi-horizontal gradient of (5.2) on constant *σ ^{O}* surfaces. This gives, after some manipulation,

For the displaced reference ocean, again *ψ ^{O}* = 1, and (5.3) pertains, so that the last factor in the second term is −

*σ*′

*Δ*

_{B}*p*; also

*N*

^{2}/

*g*≈

*ρgσ*′

*. Hence,*

_{B}This is indistinguishable from Eq. (3.6), the neutrality of neutral density, shown in Fig. 5b. One may compare it to (4.5), the neutrality of simple global orthobaric density surfaces, also shown in Fig. 5b. While the meridional dependence is shifted and reversed, the amount of variation is identical.

## 6. Summary and discussion

Density is an important variable, not merely as a tracer, but especially because of its link to motion acceleration through the buoyant force. As in situ density’s utility as a tracer is limited by its susceptibility to adiabatic compression, efforts are made to correct for this effect. This would be uncontroversial if the temperature–salinity relation of the ocean were unvarying with position. As this is not so, however, only partial solutions to the difficulty are possible. This paper has considered two attempts at such solutions: neutral density and orthobaric density. The results are judged by a number of criteria: materiality, the degree to which reversible contributions (apart from irreversible, diabatic contributions) to the substantial rate of change of the corrected density variable are minimal; neutrality, the degree to which isopleth surfaces of the variable are everywhere normal to the dianeutral vector **d** [Eq. (1.4)]; and the degree to which buoyancy forces can be regarded as acting across isopleth surfaces, embodied in the paradigm of Eq. (1.3) for isentropes when the temperature–salinity relation is unvarying. These criteria are not independent, as the resemblances of Eq. (2.10) (neutrality) to (2.9) (materiality) show.

Our approach was analytic. We constructed an idealization of the temperature–salinity structure of the upper kilometer or more of the Atlantic Ocean, in which the ocean transitions continuously with latitude from a southern *θ*–*S* relation to a northern one. A simplified equation of state was assumed, though with a representative, nonzero, thermobaric parameter, *V _{θp}*/

*V*. For the reference ocean required by Jackett and McDougall’s (1997) neutral density method, for which the Levitus climatology is usually adopted, we supposed a quiescent rest state, though with the said meridional

_{θ}*θ*–

*S*variations. On the other hand, for the reference

*θ*–

*S*relation required for simple orthobaric density, we took an intermediate relation between the northern and southern extremes. For the present ocean, as a test, we assumed a vertically displaced form (with resulting motion) of the quiescent state, and computed distributions of materiality and neutrality for each candidate variable. While this test does not exhaust the range of real-world possibilities, it throws light on the relative performances of the variables.

We compared, in Fig. 5, the material flow (materiality) across isopleths of neutral density and simple orthobaric density, and the mismatch (neutrality) of the dips of said surfaces from the neutral dip. In sections 2 and 4, we identified several contributions to each of these indices: type I, due to divergence of the local *θ*–*S* relation from the reference *θ*–*S* relation, and type II, from the meridional gradient of the reference *θ*–*S* relation. Because the meridional variation of the *θ*–*S* relation is built into the reference ocean, there is by construction no type I contribution for neutral density, though there must be a type II contribution. For simple orthobaric density, on the other hand, the reference *θ*–*S* relation is meridionally unvarying so there is no type II contribution, although, since the local *θ*–*S* relation differs from the reference relation, there must be a type I contribution. Although their spatial distributions are different, the respective types of contribution to the materiality and neutrality of the two variables are similar and comparable in magnitude. McDougall and Jackett (2005a), who reached a different conclusion, which favored neutral density, did so by overlooking the type II contributions.

A dual-hemisphere orthobaric density function was devised (section 4). It takes advantage of the close hemispheric validity of two different *θ*–*S* relations: southern and northern. This permits the local reduction of type I materiality and neutrality contributions. But a correspondence between the two distinct orthobaric density functions must be established in order to connect isopleth surfaces between hemispheres. This is achieved by supposing meridional continuity when the ocean is in the reference state. When the ocean is not in the reference state, discontinuities occur in the surfaces connected across the interhemispheric boundary (the equator). These are important because, even apart from the materiality of the surfaces in the respective hemispheres, they are sites where matter may pass through the gaps in the surfaces. This form of orthobaric density is notable because it is a model of the two-hemisphere form of patched potential density (Reid 1994) obtained in the limit when the vertical intervals (or patches, typically chosen to be 1000 dbar thick in practice) are infinitesimally small. The patched potential density has been seminal in the development of both neutral density (Jackett and McDougall 1997) and orthobaric density (de Szoeke et al. 2000; de Szoeke and Springer 2005).

The dual-hemisphere idea can be extended to form multiregional orthobaric density functions, each chosen with reference to a local *θ*–*S* relation, and joined between neighboring regions by the requirement of continuity when the ocean is in the reference state (section 5). Thus, type I materiality can be reduced virtually to nothing, while the discontinuities among regions, and the material flow across them, even though individually reduced, become more frequent. In the limit of a meridionally continuous *θ*–*S* relation, the latter material flow is seen to be an analog of the type II cross-material flow of neutral density. This analogy is quantitative. Apart from a 50% reduction factor for neutral density {which can be traced to the choice of a half-and-half reference pressure [Eq. (2.8)]}, the meridional structure of the type II materiality is identical.

Though the materiality of a variable is of obvious utility in interpreting its spatial variation, the fundamental significance of neutrality is not so clear. It is often said that the neutral plane at a specific point (the plane perpendicular to **d** at that point) contains the bundle of trajectories on which a water parcel can begin to move without change in heat or energy, or that *ρ*^{−1}(−**∇***p* + *c*^{2}**∇***ρ*) = −*c*^{2}**d** is the buoyant restoring force of adiabatic test particle excursions. We cannot substantiate such statements. Rather, when, for example, a spatially unvarying *θ*–*S* relation applies, so that neutral surfaces exist (and are coincident with isentropes), the combination of the pressure gradient force and buoyancy force with respect to isentropes may be replaced by the dynamical relation (1.3), in which the effective buoyant force is Π^{(θ)}**k**′^{(θ)} = Π^{(θ)}**∇***θ*, acting parallel to **d**. This is indeed a fundamental property of neutral surfaces when they exist. From it, familiar inferences may be drawn about potential vorticity conservation and thermal wind balances. The extension or modification of this property to the case of the spatially varying *θ*–*S* relation (so that no neutral surfaces exist) is the pertinent dynamical question. We obtained in the appendix the generalization [Eq. (A.12)], with respect to neutral density surfaces, of the substitution (1.3), and showed that there are additional terms—virtual, nonconservative forces in effect—associated with the spatial variation of the *θ*–*S* relations implicit in the reference datasets used in the construction of neutral density. These additional virtual forces are experienced by water parcels constrained to follow neutral density surfaces. A systematic study of the importance of these forces is yet lacking.

The simple orthobaric density *σ*^{⊕} is specifically constructed to exhibit under any circumstances the dynamical substitution property (1.3) that isentropes possess when the *θ*–*S* relation is unvarying. Nonsimple variants of orthobaric density will, like neutral density, exhibit extra virtual forces. For example, on dual-hemisphere orthobaric density surfaces (based on different northern and southern *θ*–*S* relations in each hemisphere of the globe), which may exhibit discontinuities at the equator, an impetus proportional to the discontinuity must be administered to a hypothetical water parcel crossing the equator to move it to the designated continuation of the surface. This is analogous to the continuous virtual forces on neutral density surfaces. Such forces must be applied to water parcels to keep them moving on continuous, extended orthobaric density surfaces.

An important question remains about the nature of orthobaric density and neutral density. It concerns the sizes of the contributions to materiality that have been the main subject of this paper—the reversible terms, principally types I and II—relative to the irreversible, adiabatic terms (denoted as type IV in sections 2 and 4 and in the appendix). The angle between simple orthobaric density surfaces and neutral planes has been estimated on long meridional hydrographic sections through the Atlantic and Pacific Oceans [see de Szoeke et al.’s (2000) Fig. 12]. This angle is to very good approximation the neutrality magnitude |**∇**_{σ⊕}*z* − **h**| calculated in section 4. De Szoeke et al. (2000) and McDougall and Jackett (2005a) reached different conclusions about the significance of the errors incurred in assuming that the diffusivity tensor’s principal axes are aligned with orthobaric density surfaces (or neutral density surfaces) rather than along neutral planes. However, a recent critical examination of these various alignment hypotheses shows that they are superfluous and makes the divergence of views irrelevant and moot (de Szoeke 2009, unpublished manuscript).

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX

#### Neutral Density

##### Materiality

The general definition of neutral density, Eq. (2.8), gives the target pressure *p ^{r}* in terms of a water sample’s values of potential temperature

*θ*and salinity

*S*at ambient pressure

*p*and horizontal location

**x**, using archived reference profiles

*θ*

_{ref}and

*S*

_{ref}at that location. Equation (2.7) with

*p*substituted then gives

^{r}*γ*. It will be convenient, though not essential, to use the following simplified equation of state:

where

Here, *ρ* is the density, *V* is the specific volume, *δ* is the SVA normalized by *V*_{0}, *θ* is the potential temperature, *S* is the salinity, *V*(*p*) = *V*(*p*,*θ*_{0},*S*_{0}), and *μ* = (*V _{pθ}*/

*V*)

_{θ}_{p=0}≈ 2 × 10

^{−4}dbar

^{−1}is the thermobaric coefficient (McDougall 1987b).

^{A1}In (A.1b), we used the parameters

*θ*

_{0}= 5°C,

*S*

_{0}= 34.63 psu,

^{A2 }

*a*

_{1}= 1.13 × 10

^{−4}°C

^{−1},

*a*

_{2}= 1.35 × 10

^{−5}°C

^{−2}, and

*b*= 0.799 × 10

^{−3}psu

^{−1}, for the first and second thermal expansion coefficients, and the haline contraction coefficient, respectively. The potential density

*ρ*with respect to a fixed reference pressure

_{r}*r*is given by

When the material derivatives of both sides of (2.8) and (2.7) are taken, and *ṗ ^{r}* ≡

*Dp*/

^{r}*Dt*=

**u**·

**∇**

*+*

_{γ}p^{r}*γ̇*∂

*is eliminated, one obtains the following relation among the material derivative*

_{γ}p^{r}*γ̇*of neutral density and

*ϑ̇*,

*ṡ*,

*ṗ*, and

**ẋ**≡

**u**:

The first equality of (A.2) defines the material velocity *e ^{γ}* across surfaces of constant neutral density

*γ*; the Roman numerals identify terms referred to in the main text; and

*ψ̂*is the buoyancy gain factor for neutral density, defined by

^{γ}The first equality in (A.4) defines *D*, while the second follows from taking the partial *γ* derivative, at constant **x** and *t*, of both sides of Eq. (2.8). Also,

Here, 𝗞 is the turbulent diffusivity tensor, *q̇*^{(p)} represents the total effect on the negative buoyancy of irreversible exchanges of heat and salt [the superscript (*p*) emphasizes that the normalized thermal expansion coefficient *V*_{0}^{−1}∂*V*/∂*ϑ* = 1 + *μ**p*, occurring in (A5), is calculated at pressure *p*], and *a*_{2} = *V*_{0}^{−1}∂^{2}*V*/∂*θ*^{2} is the second expansion coefficient, or cabbeling coefficient, defined after Eq. (A.1b). In obtaining (A.2) use was made of identities and definitions such as

Here, *ϑ*_{ref}, *s*_{ref}, *N*_{ref}, **d**_{ref} are functions of *p ^{r}* and position

**x**{with

*ρ*

_{ref}= [

*V*(

*p*,

^{r}*ϑ*

_{ref},

*s*

_{ref})]

^{−1}} in the reference ocean, to be evaluated at

*p*=

^{r}*p*

_{ref}(

*γ*,

**x**), while

*ϑ*,

*s*, and

*N*are functions of

*p*(or

*γ*),

**x**, and

*t*in the present ocean. The equalities in the second and third lines of (A.6) define the vertical and horizontal components of

**d**

_{ref}, the dianeutral vector in the reference dataset. The last equality defines the neutral dip

**h**

_{ref}, a vector whose magnitude and direction in the horizontal plane give the gradient and orientation, respectively, of neutral planes.

The derivation of (A.2) has as yet used no particular property of the reference neutral density surfaces. However, the vanishing, if possible, of the factor in parentheses of term III in (A.2) (which concerns only fields from the reference dataset) could be used to provide *p ^{r}* =

*p*

_{ref}(

*γ*,

**x**) [or its inverse

*γ*=

*γ*

_{ref}(

*p*,

^{r}**x**)]. The exact vanishing of term III is possible if and only if

**d**

_{ref}·

**∇**×

**d**

_{ref}= 0. If this condition is not fulfilled, one may settle for

*p*

_{ref}(

*γ*,

**x**), which makes term III ≈ 0 in some appropriate global minimal sense (Eden and Willebrand 1999).

Equation (A.2) is similar to an equation derived by McDougall and Jackett (2005a,b), except that in the first expression in (A.4) for the factor *D* these authors in effect assume that *N*_{ref} ≈ *N* and ∂*p ^{r}*/∂

*γ*≈ ∂

*p*/∂

*γ*. We retain the generality of (A.4) because the buoyancy frequency and spacing of neutral density surfaces may differ substantially from the reference ocean to the instantaneous present ocean. We have also retained generality in the choice of weights

*λ*and

*λ̃*in the reference pressure

*p*=

*λ̃p*+

*λp*of the potential densities in Eq. (2.8); this includes McDougall and Jackett’s (2005a,b) standard setting of

^{r}*λ*=

*λ̃*= 0.5. McDougall and Jackett’s

*ψ*factor is related to (A.3) by

^{γ}which has a limiting value even when *ϑ* → *ϑ*_{ref}.

One may compare (A.2) to (5.4) for the flow across extended orthobaric isopycnals. Term I of Eq. (A.2) is the sum of two terms weighted, respectively, by *λ̃* and *λ*, and both multiplied by *ψ̂ ^{γ}* − 1. The first term is the apparent vertical motion (∂

*+*

_{t}**u**·

**∇**)

*of neutral density surfaces*

_{γ}p*p*=

*p*(

*γ*,

**x**,

*t*) (as though they were material) in the present instantaneous dataset. The second is the apparent vertical motion

**u**·

**∇**

*of neutral density surfaces*

_{γ}p^{r}*p*=

^{r}*p*

_{ref}(

*γ*,

**x**) frozen in the reference dataset, as though the present instantaneous current

**u**=

**u**(

*γ*,

**x**,

*t*) was projected into the reference dataset on the corresponding neutral density surface (Fig. 1). This weighted sum of apparent vertical motions seems rather peculiar. Equation (5.4) contains a counterpart only of the first of these two contributions. [If the weights

*λ̃*= 1,

*λ*= 0 were chosen, the correspondence of terms between (A.2) and (5.4) would be very close.]

Most important, the second term of (5.4) for the material derivative of the extended orthobaric density corresponds to the *p* − *p ^{r}* term (II) of (A.2). For

*λ*=

*λ̃*= 0.5, using Eq. (A.3) and

*p*= −

_{z}*ρg*, the first few terms of Eq. (A.2) become

McDougall and Jackett (2005a) assert that, “Since (*p* − *p ^{r}*)

**u**·

**∇**

_{γ}ϑ_{ref}is typically of the same order as (

*ϑ*−

*ϑ*

_{ref})

**u**·

**∇**

*, all contributions to [Eq. (A.2)′] typically have the same magnitude.” Then, because*

_{γ}p^{r}*ϑ*−

*ϑ*

_{ref}∝

*ψ*− 1 ≈ 0, it would follow that all shown terms, types I and II, are negligible. However, in the main text we contest the quoted statement.

^{γ}##### Neutrality

By taking the gradient along neutral density surfaces **∇*** _{γ}* of both sides of Eq. (2.8), and making use of the definitions and identities (A.4) and (A.6), one arrives at

(The Roman numerals refer to terms discussed in the main text.) The two quantities in the parentheses on the left side and in term III on the right of (A.9) are the neutralities in the present and reference datasets, that is, the respective differences between the neutral density surface dip and the neutral plane dip. Even if the reference neutrality on the right of (A.9) were zero, the other two terms, particularly the second, proportional to *μ*(*p*−*p ^{r}*), may significantly bias the neutrality in the present dataset.

##### Generalized Montgomery function

or

where

One may regard Eq. (A.10) as a metaequation of state, furnishing density as a function of metastate variables *γ* and **x**, in addition to *p* and *ϑ*. Defining a generalized heat content

and a generalized Montgomery function *M*^{(γ)} = *H*^{(γ)} + *gz*, one can show that

where **k**′^{(γ)} = (**k** − **∇*** _{γ}z*)/∂

*=*

_{γ}z**∇**

*γ*, and

Here Π^{(γ)}, Θ, and **X** are generalized metathermodynamic forces given by the respective gradients of *H*^{(γ)}. Equation (A.12) is the generalization of the dynamical substitution property (1.3). To quantitatively consider distributions of hydrographic variables on neutral density surfaces and the motion fields that might have brought them about, one must take account of the extra forces that occur in the momentum balances with respect to such surfaces. For this the generalized substitution (A.12) is relevant and indispensable.

## Footnotes

* Current affiliation: Earth and Space Research, Seattle, Washington.

*Corresponding author address:* Roland A. de Szoeke, 104A COAS Administration Bldg., Oregon State University, Corvallis, OR 97331. Email: szoeke@coas.oregonstate.edu

^{1}

It has been traditional in physical oceanography to denote as conservative what we call here material. Because the former has recently been used in a rather different sense (McDougall 2003), we have adopted the latter term.

^{2}

Notation conventions: **∇***p* = **∇*** _{z}p* +

**k**∂

*and*

_{z}p**∇**

*=*

_{z}p**i**∂

_{x|z}

*p*+

**j**∂

_{y|z}

*p*, where ∂

_{x|z}

*p*=

*p*

_{x|z}, ∂

_{y|z}

*p*=

*p*

_{y|z}, etc., are partial derivatives at fixed

*z*; ∂

*=*

_{z}p*p*, etc., are partial derivatives with respect to

_{z}*z*;

**i**,

**j**,

**k**are orthogonal unit basis vectors, horizontal and vertical;

**∇**

*=*

_{V}z**i**∂

_{x|V}

*z*+

**j**∂

_{y|V}

*z*is the gradient of isopycnal surfaces; ∂

_{x|V}

*z*=

*z*

_{x|V}, ∂

_{y|V}

*z*=

*z*

_{y|V}, etc., are partial derivatives at constant

*V*; ∂

*=*

_{V}z*z*, etc., are partial derivatives with respect to

_{V}*V*; and

**∇**

_{V}M^{(V)}=

**i**∂

_{x|V}

*M*

^{(V)}+

**j**∂

_{y|V}

*M*

^{(V)}is the two-dimensional vector gradient at constant

*V*of the Montgomery function

*M*

^{(V)}. Similar usages with independent variables other than

*V*will be employed.

^{3}

Neutral density and its isopleth surfaces are well defined, even when neutral surfaces do not exist (i.e., nonzero helicity). We will be careful to refer to neutral density surfaces, not neutral surfaces.

^{4}

Surfaces may be labeled alternatively by any Γ(*γ*). It is convenient to resolve this ambiguity by setting *γ* to the potential density at a central location (e.g., Jackett and McDougall 1997).

^{5}

We prefer to use pressure *p* rather than height *z* as a vertical coordinate.

^{6}

The material derivative operator is given by any of the following: *D*/*Dt* = ∂_{t|z} + **u** · **∇**_{z} + *w*∂_{z} = ∂_{t|p} + **u** · **∇**_{p} + *ṗ*∂_{p} = ∂_{t|γ} + **u** · **∇**_{γ} + *γ̇*∂_{γ}, etc., where *w* = *ż* = *Dz*/*Dt*, *ṗ* = *Dp*/*Dt*, *γ̇* = *Dγ*/*Dt*, etc.; ∂_{t|z}, etc., are partial time derivatives with *z*, etc., held constant. This invariance to reference frame is demonstrated by de Szoeke and Samelson (2002).

^{7}

McDougall and Jackett (2005a) define *ψ ^{γ}* − 1 =

*λ*(

*ψ̂*− 1), which is convenient when

^{γ}*λ*=

*λ̃*= 0.5 is used.

^{9}

More systematically, de Szoeke et al. (2000) defined a global orthobaric density function on the basis of a volumetrically averaged compressibility function—equivalent to a volumetrically averaged *θ*–*S* relation.

^{A1} More generally, what follows may be readily extended to cover any equation of state of the form *δ* = *P*(*p*)*ϑ*(*θ*, *S*) − *s*(*θ*, *S*), where *θ* and *S* are the true potential temperature and salinity.

^{A2} Here, psu refers to United Nations Educational, Scientific and Cultural Organization (UNESCO) Practical Salinity Scale of 1978 (PSS78) salinity units (UNESCO 1981).