Neutral density

The neutral density ( formula_1 ) or empirical neutral density is a density variable used in oceanography, introduced in 1997 by David R. Jackett and Trevor McDougall.
It is 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 formula_1 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 ("dianuetral 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. 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 demonstrated that the neutral tangent plane, and hence also neutral surfaces, are normal to the dianeutral vector
where formula_4 is the salinity, formula_5 is the potential temperature, formula_6 the thermal expansion coefficient and formula_7 the saline concentration coefficient.
Thus, neutral surfaces are defined as surfaces everywhere perpendicular to formula_8.
The contribution to density caused by gradients of formula_4 and formula_10 within the surface exactly compensates. That is, with formula_11 the 2D gradient within the neutral surface,
If such a neutral surface exists, the neutral helicity formula_13 (related in form to hydrodynamical helicity) must be zero everywhere on that surface, a condition arising from non-linearity of the equation of state. 
A continuum of such neutral surfaces could be usefully represented as isosurfaces of a 3D scalar field formula_14 that satisfies
if the residual formula_16. Here, formula_17 is an integrating scalar factor that is function of space.
A necessary condition for the existence of formula_14 with formula_19 is that formula_20 everywhere in the ocean. However, islands complicate the topology such that this is not a sufficient condition. 
In the real ocean, the neutral helicity formula_21 is generally small but not identically zero. Therefore, it is impossible to create analytically a well-defined neutral surfaces, nor a 3D neutral density variable such as formula_22. 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 formula_23. Similarly, it is only possible to define formula_14 satisfying () with formula_25. Numerical techniques can be used to solve the coupled system of first-order partial differential equations () while minimizing some norm of formula_26. 
Jackett and McDougall provided such a formula_1 having small formula_28, and demonstrated that the inaccuracy due to the non-exact neutrality (formula_29) is below the present instrumentation error in density. Neutral density surfaces stay within a few tens meters of an ideal neutral surface anywhere in the world.
Given how formula_1 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 and ).
Spatial dependence.
Neutral density is a function of latitude and longitude. This spatial dependence is a fundamental property of neutral surfaces. From (), the gradients of formula_4 and formula_10 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 formula_10 per formula_10 value (or, equivalently expressed by formula_4). Thus, the connectedness of level sets of formula_10 on a neutral surface is a vital topological consideration. These regions are precisely those regions associated with the edges of the Reeb graph of formula_10 on the surface, as shown by Stanley.
Given this spatial dependence, calculating neutral density requires knowledge of the spatial distribution of temperature and salinity in the ocean. Therefore, the definition of formula_1 has to be linked with a global hydrographic dataset, based on the climatology of the world’s ocean (see World Ocean Atlas and ).
In this way, the solution of () provides values of formula_1 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 () can be solved over a more limited set of data.
Algorithm for the computation of neutral surfaces using formula_1.
Jackett and McDougall constructed the variable formula_1 using the data in the “Levitus dataset”.
As this dataset consists of measurements of S and T at 33 standard depth levels at a 1° resolution, the solution of () 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 () 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 () are neutral surfaces along which formula_1 is constant) and the finite differences method in the remaining 15%.
The output of these calculations is a global dataset labeled with values of formula_1.
The field of formula_1 values resulting from the solution of the differential system () satisfies () an order of magnitude better (on average) than the present instrumentation error in density.
The labeled dataset is then used to assign formula_1 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 formula_1.
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.
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 formula_1 values to any arbitrary hydrographic data.
Another function provided in the code, given a vertical profile of labeled data and formula_1 surfaces, finds the positions of the specified formula_1 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 formula_1 and the previous commonly used methods to obtain discretely referenced neutral surfaces (see for example Reid (1994), 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) and an easier and computationally less expensive algorithm to form neutral surfaces. 
A neutral surface defined using formula_1 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. If potential density surfaces are used, the difference can be hundreds of meters, a far larger error.