The neutral density ( ${\displaystyle \gamma ^{n}\,}$ ) or empirical neutral density is a density variable used in oceanography, introduced in 1997 by David R. Jackett and Trevor McDougall.^[1] It is a function of the three state variables (salinity, temperature, and pressure) and the geographical location (longitude and latitude). It has the typical units of density (M/V). Isosurfaces of ${\displaystyle \gamma ^{n}\,}$ form “neutral density surfaces”, which are closely aligned with the "neutral tangent plane". It is widely believed, although this has yet to be rigorously proven, that the flow in the deep ocean is almost entirely aligned with the neutral tangent plane, and strong lateral mixing occurs along this plane ("epineutral mixing") vs weak mixing across this plane ("dianeutral mixing"). These surfaces are widely used in water mass analyses. Neutral density is a density variable that depends on the particular state of the ocean, and hence is also a function of time, though this is often ignored. In practice, its construction from a given hydrographic dataset is achieved by means of a computational code (available for Matlab and Fortran), that contains the computational algorithm developed by Jackett and McDougall. Use of this code is currently restricted to the present day ocean.

Mathematical expression

The neutral tangent plane is the plane along which a given water parcel can move infinitesimally while remaining neutrally buoyant with its immediate environment.^[1] This is well-defined at every point in the ocean. A neutral surface is a surface that is everywhere parallel to the neutral tangent plane. McDougall^[2] demonstrated that the neutral tangent plane, and hence also neutral surfaces, are normal to the dianeutral vector

${\displaystyle \mathbf {N} =\rho (\beta \nabla S-\alpha \nabla \theta ),}$

where ${\displaystyle S}$ is the salinity, ${\displaystyle \theta \,}$ is the potential temperature, ${\displaystyle \alpha \,}$ the thermal expansion coefficient and ${\displaystyle \beta \,}$ the saline concentration coefficient. Thus, neutral surfaces are defined as surfaces everywhere perpendicular to ${\displaystyle \mathbf {N} }$. The contribution to density caused by gradients of ${\displaystyle S}$ and ${\displaystyle \theta }$ within the surface exactly compensates. That is, with ${\displaystyle \nabla _{n))$ the 2D gradient within the neutral surface,

${\displaystyle \beta \nabla _{n}S=\alpha \nabla _{n}\theta .}$ (1)

If such a neutral surface exists, the neutral helicity ${\displaystyle H=\mathbf {N} \cdot \nabla \times \mathbf {N} }$ (related in form to hydrodynamical helicity) must be zero everywhere on that surface, a condition arising from non-linearity of the equation of state.^[3] A continuum of such neutral surfaces could be usefully represented as isosurfaces of a 3D scalar field ${\displaystyle \gamma ^{n))$ that satisfies^[1]

${\displaystyle \nabla \gamma ^{n}\ =b\mathbf {N} +{\cal {R)),}$ (2)

if the residual ${\displaystyle {\cal {R))=0}$. Here, ${\displaystyle b}$ is an integrating scalar factor that is function of space.

A necessary condition for the existence of ${\displaystyle \gamma ^{n))$ with ${\displaystyle {\cal {R))=0}$ is that ${\displaystyle H=0}$ everywhere in the ocean.^[1] However, islands complicate the topology such that this is not a sufficient condition.^[4]

In the real ocean, the neutral helicity ${\displaystyle H}$ is generally small but not identically zero.^[5] Therefore, it is impossible to create analytically a well-defined neutral surfaces, nor a 3D neutral density variable such as ${\displaystyle \gamma ^{n))$.^[6] There will always be flow through any well-defined surface caused by neutral helicity.

Therefore, it is only possible to obtain approximately neutral surfaces, which are everywhere _approximately_ perpendicular to ${\displaystyle \mathbf {N} }$. Similarly, it is only possible to define ${\displaystyle \gamma ^{n))$ satisfying (2) with ${\displaystyle {\cal {R))\neq 0}$. Numerical techniques can be used to solve the coupled system of first-order partial differential equations (2) while minimizing some norm of ${\displaystyle {\cal {R))}$.

Jackett and McDougall^[1] provided such a ${\displaystyle \gamma ^{n}\,}$ having small ${\displaystyle {\cal {R))}$, and demonstrated that the inaccuracy due to the non-exact neutrality (${\displaystyle {\cal {R))\neq 0}$) is below the present instrumentation error in density.^[7] Neutral density surfaces stay within a few tens meters of an ideal neutral surface anywhere in the world.^[8]

Given how ${\displaystyle \gamma ^{n}\,}$ has been defined, neutral density surfaces can be considered the continuous analog of the commonly used potential density surfaces, which are defined over various discrete values of pressures (see for example ^[9] and ^[10]).

Spatial dependence

Neutral density is a function of latitude and longitude. This spatial dependence is a fundamental property of neutral surfaces. From (1), the gradients of ${\displaystyle S}$ and ${\displaystyle \theta }$ within a neutral surface are aligned, hence their contours are aligned, hence there is a functional relationship between these variables on the neutral surface. However, this function is multivalued. It is only single-valued within regions where there is at most one contour of ${\displaystyle \theta }$ per ${\displaystyle \theta }$ value (or, equivalently expressed by ${\displaystyle S}$). Thus, the connectedness of level sets of ${\displaystyle \theta }$ on a neutral surface is a vital topological consideration. These regions are precisely those regions associated with the edges of the Reeb graph of ${\displaystyle \theta }$ on the surface, as shown by Stanley.^[4]

Given this spatial dependence, calculating neutral density requires knowledge of the spatial distribution of temperature and salinity in the ocean. Therefore, the definition of ${\displaystyle \gamma ^{n}\,}$ has to be linked with a global hydrographic dataset, based on the climatology of the world's ocean (see World Ocean Atlas and ^[11]). In this way, the solution of (2) provides values of ${\displaystyle \gamma ^{n}\,}$ for a referenced global dataset. The solution of the system for a high resolution dataset would be computationally very expensive. In this case, the original dataset can be sub-sampled and (2) can be solved over a more limited set of data.

Algorithm for the computation of neutral surfaces using neutral density

Jackett and McDougall constructed the variable ${\displaystyle \gamma ^{n}\,}$ using the data in the "Levitus dataset".^[12] As this dataset consists of measurements of S and T at 33 standard depth levels at a 1° resolution, the solution of (2) for such a large dataset would be computationally very expensive. Therefore, they sub-sampled the data of the original dataset onto a 4°x4° grid and solved (2) on the nodes of this grid. The authors suggested to solve this system by using a combination of the method of characteristics in nearly 85% of the ocean (the characteristic surfaces of (2) are neutral surfaces along which ${\displaystyle \gamma ^{n}\,}$ is constant) and the finite differences method in the remaining 15%. The output of these calculations is a global dataset labeled with values of ${\displaystyle \gamma ^{n}\,}$. The field of ${\displaystyle \gamma ^{n}\,}$ values resulting from the solution of the differential system (2) satisfies (2) an order of magnitude better (on average) than the present instrumentation error in density.^[13]

The labeled dataset is then used to assign ${\displaystyle \gamma ^{n}\,}$ values to any arbitrary hydrographic data at new locations, where values are measured as a function of depth by interpolation to the four closest points in the Levitus atlas.

Practical computation of neutral density

The formation of neutral density surfaces from a given hydrographic observation requires only a call to a computational code that contains the algorithm developed by Jackett and McDougall.^[14]

The Neutral Density code comes as a package of Matlab or as a Fortran routine. It enables the user to fit neutral density surfaces to arbitrary hydrographic data and just 2 MBytes of storage are required to obtain an accurately pre-labelled world ocean.

Then, the code permits to interpolate the labeled data in terms of spatial location and hydrography. By taking a weighted average of the four closest casts from the labeled data set, it enables to assign ${\displaystyle \gamma ^{n}\,}$ values to any arbitrary hydrographic data.

Another function provided in the code, given a vertical profile of labeled data and ${\displaystyle \gamma ^{n}\,}$ surfaces, finds the positions of the specified ${\displaystyle \gamma ^{n}\,}$ surfaces within the water column, together with error bars.

Advantages of using the neutral density variable

Comparisons between the approximated neutral surfaces obtained by using the variable ${\displaystyle \gamma ^{n}\,}$ and the previous commonly used methods to obtain discretely referenced neutral surfaces (see for example Reid (1994),^[10] that proposed to approximate neutral surfaces by a linked sequence of potential density surfaces referred to a discrete set of reference pressures) have shown an improvement of accuracy (by a factor of about 5) ^[15] and an easier and computationally less expensive algorithm to form neutral surfaces. A neutral surface defined using ${\displaystyle \gamma ^{n}\,}$ differs only slightly from an ideal neutral surface. In fact, if a parcel moves around a gyre on the neutral surface and returns to its starting location, its depth at the end will differ by around 10m from the depth at the start.^[8] If potential density surfaces are used, the difference can be hundreds of meters, a far larger error.^[8]