A new method for forming approximately neutral surfaces

A new method for forming approximately neutral surfaces A. Klocker, T. J. McDougall, and D. R. Jackett Centre for Australian Climate and Weather Research, Castray Esplanade, TAS 7000, Australia Antarctic Climate and Ecosystems Cooperative Research Center, University of Tasmania, Private Bag 80, TAS 7001, Australia Institute of Antarctic and Southern Ocean Studies, University of Tasmania, Private Bag 77, TAS 7001, Australia Received: 11 July 2008 – Accepted: 14 July 2008 – Published: 29 August 2008 Correspondence to: A. Klocker (andreas.klocker@csiro.au) Published by Copernicus Publications on behalf of the European Geosciences Union.


Introduction
Transport in the ocean does not occur along surfaces of constant in situ density and several approaches have been used to find a density variable whose isosurfaces accurately describe the direction along which flow in the ocean occurs.Using inappropriate density surfaces leads to a fictitious diapycnal diffusivity, D f , sometimes orders of magnitude larger than the measured diapycnal diffusivity in the ocean.D f is an error resulting from mixing along a well-defined surface instead of along neutral tangent planes.This fictitious diapycnal diffusivity does not represent a real physical process.
Correspondence to: A. Klocker (andreas.klocker@csiro.au)Scalar properties in the ocean get stirred (and subsequently mixed) efficiently by mesoscale eddies and two-dimensional turbulence along neutral tangent planes (McDougall, 1987).These are defined such that when water parcels are moved small distances along these planes, they experience no buoyant restoring forces.It is impossible to link these neutral tangent planes to form a surface, therefore "neutral surfaces" will always be mathematically ill-defined (McDougall, 1987;McDougall and Jackett, 1988).If we were to follow a neutral trajectory around an ocean basin (linking up neutral tangent planes) and arrive back at the initial latitude/longitude one normally arrives at a different depth than where one started.This shows that the definition of a neutral surface is pathdependent, an effect caused by the nonlinearity of the equation of state of seawater (because the ratio α /β is a function of pressure; see Appendix A for a more detailed explanation).Therefore it is not possible to find a "perfect" surface to describe flow in the ocean.There will always be errors associated with density surfaces due to path-dependency -but how large is this unavoidable error?Efforts to construct density variables minimizing D f include approximately neutral surfaces (Jackett and Mc-Dougall, 1997;Jackett et al., 2009) and orthobaric density surfaces (de Szoeke et al., 2000).These algorithms label a three-dimensional hydrography with a density variable.We can then find surfaces in this hydrography on which the density variable is constant and use this surface for inverse techniques or for plotting variables such as temperature, salinity and nutrients to understand the evolution of water masses.Compared to these density-labelling algorithms the technique described in this work takes one density surface -which can be a surface of constant potential density, neutral density or any other density variable -and improves it to ensure it is as close to the neutral tangent planes as possible thus minimizing the fictitious diapycnal diffusivity.This algorithm is ideal for creating optimized approximately neutral Published by Copernicus Publications on behalf of the European Geosciences Union.
A. Klocker et al.: Optimized approximately neutral surfaces surfaces to use as water mass density boundaries in inverse models.The improvement due to these optimized approximately neutral surfaces might not be significant in large box inverse models of non-synoptic hydrographic sections compared to the other assumptions made (i.e.steady state, etc.) but we expect that these surfaces will significantly decrease the error of inverse models using synoptic sections for process studies that particularly target the determination of mixing.

Basic properties of density surfaces
Many different density surfaces have been used in the past for inverse techniques, layered ocean models or other applications describing ocean circulation along isopycnals.These different density surfaces all differ in the extent to which they achieve the three desirable but mutually inconsistent properties (McDougall and Jackett, 2005b): -being as neutral as possible -being as quasi-material as possible -possessing a geostrophic streamfunction (commonly called a Montgomery potential), where quasi-material means that flow through a surface only arises due to mixing processes.
In this work we will mainly focus on the first point, comparing how "neutral" different density variables are."Neutral" here describes the direction along which a parcel can travel without experiencing buoyant restoring forces.The γ n -variable (Jackett and McDougall, 1997) and the γ ivariable (Jackett et al., 2009) for example were constructed to produce a surface which is as neutral as possible by minimizing the slope difference between these respective surfaces and the neutral tangent planes, but ignoring the last two points mentioned above.Eden and Willebrand (1999) took a different approach and tried to construct a density variable which is a compromise between neutrality and two other properties, (a) the horizontal gradient of the neutral density should agree with the gradient of in situ density and (b) the vertical gradient of the neutral density should be proportional to the static stablility of the water column.
These requirements are quite different to the properties used by McDougall and Jackett (2005b).We note that the integrating factor b (McDougall and Jackett, 1988), defined by varies in the ocean whereas the extra requirements of Eden and Willebrand (1999) would only be strictly true if the integrating factor b were equal to one everywhere in the ocean.ρ in this equation is in situ density, β is the saline contraction coefficient, α the thermal expansion coefficient and S z and z are the vertical gradients of salinity and conservative temperature.
McDougall (1988) shows that b is given by b where N 2 is the buoyancy frequency and ∇ a is the gradient along an approximately neutral surface.T b is the thermobaric parameter given by This equation was actually derived for spatial gradients along a neutral tangent plane and here it is written in terms of gradients in an approximately neutral surface.It was also derived ignoring the dependence of the saline contraction coefficient on pressure (in comparison to α p ).For both these reasons we use an approximately equal sign in Eq. ( 2).Choosing the appropriate density variable will always depend on ones application -a surface which satisfies all three properties does not exist due to the nature of the equation of state.It is therefore very important to know the advantages of each density variable and the errors associated with them.One density variable might do a good job for one application but introduce substantial errors for another.
To quantify the quality (in the sense of being close to neutral) of a density surface we use the fictitious diapycnal diffusivity of density caused by mixing laterally along a density surface with a slope different to that of the neutral tangent plane.D f is given by where K is a lateral diffusivity (taken to be 1000 m 2 s −1 in the following calculations) and s is the slope difference between the density surface used and the neutral tangent plane, where ∇ n is the gradient along a neutral tangent plane and ∇ a is the gradient along any approximate surface (whether it be a potential density surface, an approximately neutral surface or any other surface).The fictitious diapycnal diffusivity described here is the same as D fictitious in McDougall and Jackett, 2005b.A derivation of D f can be found in Appendix B.
It has been shown that the mean diapycnal diffusivity in the ocean is roughly 10 −5 m 2 s −1 , even though it can be larger above rough topography (Polzin et al., 1997).If D f for a specific density surface is comparable or larger than Ocean Sci., 5, 155-172, 2009 www.ocean-sci.net/5/155/2009/ as is illustrated in Fig. 1, where γ a is the variable which is constant in the approximately neutral surface.The property of the ocean's hydrography which stops us from forming mathematically well-defined neutral surfaces, with N2 being the buoyancy frequency and g the gravitational acceleration.In the neutral tangent plane ǫ = 0.
An important relationship in the neutral framework is that between neutral helicity in an approximately neutral surface and the two-dimensional curl of ǫ, ∇ a × ǫ.According to theory (Eqns. (38) and (39) in McDougall and Jackett (1988)) they should be related as can be seen from the following equation Here δz is the depth change of a neutral trajectory after completing a closed loop around the ocean (i.e. the neutral trajectory finishes at the same horizontal position as it Fig. 1.This figure shows the differences between e, e hel , e a and w.The lateral velocity (ui+vj ) is directed horizontally.The three surfaces shown are the approximately neutral surface (γ a ), the neutral tangent plane (ntp) and the top-most is the lateral velocity plus w, where w includes all components leading to a flow which differs from a purely horizontal flow (the tilt of an approximately neutral surface, mixing effects and a diapycnal velocity caused by the illdefined nature of neutral surfaces, e hel ). 10 −5 m 2 s −1 over a significant area, then using this surface to describe the flow in the ocean would introduce significant mixing that is purely due to the error of the definition of the density surface used.
When describing ocean flow the terms "isopycnal" and "diapycnal" are used to describe flow along and through "density" surfaces, respectively.But since it is impossible to construct a mathematically well-defined neutral surface due to the effects of the nonlinear equation of state, it is only possible to define an approximately neutral surface.Therefore to properly define diapycnal transport we have to distinguish between flow across a mathematically well-defined approximately neutral surface and the actual isopycnal/diapycnal transport.In the latter "isopycnal" means along a neutral helix (the trajectory we would get if we connect neutral tangent planes following fluid flow) and "diapycnal" means across this neutral helix.
The vertical velocity e a through an approximately neutral surface, γ a , can be written as where e is the diapycnal transport due to cabbeling, thermobaricity, double diffusion and small-scale turbulent mixing, and e hel is the vertical velocity through the approximately neutral surface due to the helical shape of neutral trajectories (see Klocker and McDougall (2009) for an estimate of e hel ).The diapycnal velocity e hel transports mass, salinity, conservative temperature and all other tracers.This diapycnal transport, e hel , exists without requiring the dissipation of kinetic energy.It can be written as where V is the horizontal velocity (u, v).we can use a direct inversion, Alternatively we could solve Eqn.( 14) using an iterative technique, e.g., the LSQR algorithm of Paige and Saunders (1982), as implemented in Matlab (2007).For the larger data sets the iterative technique is the computationally more efficient approach.
We now have a Φ ′ -field which we need to convert into a depth change, δz, to find the depth of the optimized approximately neutral surface.From McDougall and Jackett (1988) we know that We thus have to calculate N 2 on the surface to find the depth of the new optimized approximately neutral surface 2 .
Due to the algorithm working on an horizontally extensive two-dimensional surface the first guess of the depth change, δz, will not be the final solution.We thus linearly interpolate S and Θ onto the new surface and calculate new lateral density gradient errors.We then treat these new density gradient errors as we did the ǫ init -field before to get a more accurate optimized approximately neutral surface.If we repeat these steps often enough |ǫ| 2 will converge.Once it has converged the surface will be as close to neutral as possible, with the residual fictitious diapycnal diffusivity being due only to the path-dependency caused by neutral helicity.We will call this surface the ω-surface.
On all density surfaces we will have regions where the surface outcrops or hits the bottom topography.Due to this we will end up with several regions on a density surface which do not communicate with each other.A typical example is a marginal sea with narrow connections to the open ocean.In the algorithm described above we deal with this problem of independent regions by writing a set of equations as in Eqn.( 14) for each seperate region.Similar to before we constrain the average pertubation density of each region, Φ ′ , to be zero.
When optimizing approximately neutral surfaces with the method above we sometimes get a result where |ǫ| 2 does not converge.This is because of the algorithm overestimating the depth change, δz, due to the algorithm not knowing about the stratification above and below the surface optimized.This can then lead to a growing |ǫ| 2 due to the algorithm trying to overcorrect at these casts.If this happens we have to dampen the depth change; this means we only use a certain percentage of the depth change estimated by the algorithm to calculate the optimized surface.Another way of minimizing the possibility of this problem is to discard the data in the mixed layer -a region in which other processes than neutral physics Fig. 2. The blue surface shows the initial approximately neutral surface (γ a -surface) on which the density pertubation field, , is calculated.This density pertubation field is converted into a depth change, δz, which is then applied to the initial surface to get the new apaproximately neutral surface (the blue surface).
The diapycnal transport e a can also be written in terms of the material derivative of γ a , as is illustrated in Fig. 1, where γ a is the variable which is constant in the approximately neutral surface.
The property of the ocean's hydrography which stops us from forming mathematically well-defined neutral surfaces, neutral helicity, can be written as (McDougall and Jackett, 1988): From this equation we can see that neutral helicity is a consequence of the thermobaric parameter (Eq.3), therefore a consequence of the equation of state of seawater being nonlinear in the sense that the ratio of the thermal expansion coefficient to the saline contraction coefficient is a function of pressure.
The first part of Eq. ( 9) means that for neutral helicity to be zero the line of intersection of the S and planes, ∇S×∇ , must lie in the isobaric surface, the second part requires the epineutral gradients of pressure and temperature to be parallel.Both of these requirements are close to being met in the real ocean, but the amount by which neutral helicity is non-zero may be important for some effects.
To improve existing surfaces we construct an algorithm with the aim of reducing the residual fictitious diapycnal diffusivity so that it is only due to neutral helicity and not due to any other effects.
We take one of the existing density surfaces as the initial condition and use a least-squares approach to minimize the www.ocean-sci.net/5/155/2009/Ocean Sci., 5, 155-172, 2009 ′ , is calculated.This density pertubation field is converted into a depth change, δz, which is then applied o get the new apaproximately neutral surface (the blue surface).we can use a direct inversion, Alternatively we could solve Eqn.( 14) using an iterative echnique, e.g., the LSQR algorithm of Paige and Saunders 1982), as implemented in Matlab (2007).For the larger data ets the iterative technique is the computationally more effiient approach.We now have a Φ ′ -field which we need to convert into a epth change, δz, to find the depth of the optimized approxiately neutral surface.From McDougall and Jackett (1988) e know that We thus have to calculate N 2 on th depth of the new optimized approximat Due to the algorithm working on an h two-dimensional surface the first guess δz, will not be the final solution.We thu S and Θ onto the new surface and calcu sity gradient errors.We then treat these errors as we did the ǫ init -field before to optimized approximately neutral surfac steps often enough |ǫ| 2 will converge.O the surface will be as close to neutral residual fictitious diapycnal diffusivity path-dependency caused by neutral heli surface the ω-surface.
On all density surfaces we will have r face outcrops or hits the bottom topogr will end up with several regions on a d do not communicate with each other.a marginal sea with narrow connection In the algorithm described above we de of independent regions by writing a s Eqn.( 14) for each seperate region.Simi strain the average pertubation density o be zero.
When optimizing approximately neu method above we sometimes get a resu converge.This is because of the algorith depth change, δz, due to the algorithm n stratification above and below the surf can then lead to a growing |ǫ| 2 due to th overcorrect at these casts.If this happen the depth change; this means we only u age of the depth change estimated by th late the optimized surface.Another w possibility of this problem is to discard layer -a region in which other processe 2 We use (N 2 + 3 * 10 −6 ) instead of N 2 rithm is stable when N 2 is close to zero.area integral of 2 , where is similar to the slope error s but is also dependent on vertical stratification: with N 2 being the buoyancy frequency and g the gravitational acceleration.In the neutral tangent plane =0.
An important relationship in the neutral framework is that between neutral helicity in an approximately neutral surface and the two-dimensional curl of , ∇ a × .According to theory (Eqs. 38 and 39 in McDougall and Jackett, 1988) they should be related as can be seen from the following equation Here δz is the depth change of a neutral trajectory after completing a closed loop around the ocean (i.e. the neutral trajectory finishes at the same horizontal position as it started from but at a different depth).ρ l is the locally referenced potential density.The step from −δ z N 2 g −1 to − A • dl has been derived in McDougall and Jackett (1988), and A ∇ a × •k dxdy follows from A • dl using Stokes' theorem (see Appendix C for a proof of Stokes' theorem for the two-dimensional curl).−∇ a × •k is not exactly equal to T b ∇ a p×∇ a •k because the small term β p β ∇ a p× has been ignored ( To the extent that ∇ a p and ∇ a are good approximations of ∇ n p and ∇ n , Eq. ( 9) demonstrates the approximate equivalence of T b ∇ a p×∇ a k and gN −2 H n in Eq. ( 11).
By considering a variety of areas, A, the equality of the various area integrals implies that the integrands Since neutral helicity is a property of the ocean's hydrography and we also know that −∇ a × •k is effectively equal to neutral helicity, we therefore know that −∇ a × •k is also set by the ocean's hydrography.
To check this relationship between gN −2 H n and −∇ a × • k we choose an approximately neutral surface in the North Atlantic which is close to the depth of the Mediterranean outflow (γ n =27.25 kg m −3 ).We choose this depth because one would think that this warm and salty water would cause the ocean to have increased values of neutral helicity in this region due to high temperature gradients crossing pressure gradients (see Eq. 11), making it an interesting region for our calculations of the mean diapycnal advection caused by these larger values of neutral helicity.
The data we use here and in all the following examples are model output from a standard MOM4 run with a resolution of 1 • ×2 • .The only change to the standard run is the use of conservative temperature (McDougall, 2003), , instead of potential temperature.This change is not relevant to the results.

Improvement of approximate density surfaces
Our aim is to minimize the difference between the neutral tangent planes and the approximately neutral surfaces, that is, essentially to minimize the area integral of the density gradient error .Since the curl of , ∇ a × , is given by the hydrography, we choose to minimize by adding a pertubation density field, (where =lnρ l , ρ l being the locally referenced potential density), so that is minimized while ∇ a × is unaffected by the presence of .In this way a new density surface can be formed by taking into account the pertubation density, .As described below, the new height of the density surface is adjusted by converting the density change from the density pertubation field into a depth change (using Eq. 16).This can be seen in Fig. 2 in which the blue surface is the initial approximately neutral surface (γ a -surface) on which the density pertubation field, , is calculated and the Ocean Sci., 5, 155-172, 2009 www.ocean-sci.net/5/155/2009/green surface is the new surface after the density pertubation field in the initial approximately neutral surface is imposed, where is the smallest possible density gradient error (i.e. the residual density gradient error is due to neutral helicity) and init is the initial density gradient error field.A more detailed description of the theory behind this algorithm and numerical testcases can be found in Appendix D. Now we apply the idea of minimizing , without changing its curl, ∇ a × , to construct an algorithm which optimizes existing density surfaces to be as neutral as possible with the residual error only being due to neutral helicity.As an initial condition we can use any density variable that labels a threedimensional data set.We then choose a surface on which this density variable is constant and linearly interpolate S and onto that surface.With these variables we can then calculate the density gradient error =β ∇ a S − α ∇ a1 From every grid point we want to calculate an x-component, init ew , and a y-component, init ns , of the initial density gradient error init , which we will then use as initial conditions in the algorithm.From a numerical perspective this will look like where the thermal expansion coefficient α and the saline contraction coefficient β are averaged onto the points in between the tracer grid points (the green points in Fig. 3; the red point are the tracer points).
We now construct a matrix A with the number of rows being the number of equations and the number of columns being the number of grid points.This matrix is a sparse matrix; for the init ew equations it will have a "1" for the eastern grid point and a "−1" for the western grid point -all the other entries are "0" in each row.The same is true for the init nsequations.We also constrain the average pertubation density, , to be zero.This would show up in the matrix A as a row filled with ones and in the vector init as zero.Now we have a sparse matrix A, a vector init (which has as many entries as the matrix A has rows) and we want to find the pertubation we can use a direct inversion,  One can see that D f decreases by a few orders of magnitude when using the ω-surface compared to the γ n -surface.The large improvement is possible because the model output we are using had water masses that deviated significantly from observed ocean properties.When applied to atlas data the fictitious diapycnal diffusivity in an ω-surface is perhaps just one to two orders of magnitude less than in a γ n -surface.
The improvement made by the algorithm can be seen by plotting γ n on an ω-surface (Fig. 8).γ n values have a range from 27.61 to 27.65 on the ω-surface which is a substantial density change.
Another way of seeing the improvement is by plotting  Alternatively we could solve Eq. ( 14) using an iterative technique, e.g., the LSQR algorithm of Paige and Saunders (1982), as implemented in Matlab (2007).For the larger data sets the iterative technique is the computationally more efficient approach.
We now have a -field which we need to convert into a depth change, δz, to find the depth of the optimized approximately neutral surface.From McDougall and Jackett (1988) we know that We thus have to calculate N2 on the surface to find the depth of the new optimized approximately neutral surface 2 .
Due to the algorithm working on an horizontally extensive two-dimensional surface the first guess of the depth change, Even though the change between our initial condition and the ω-surface are quite large in terms of the fictitious diapycnal diffusivity, the correlation between gN −2 H n and −∇ a ×ǫ•k or the variations of γ n on the ω-surface or the actual changes of temperature and pressure between the initial condition and the ω-surface are reasonably small.

How 'neutral' are existing density variables?
To show the differences between different density variables we use a density surface in the North Atlantic with an average depth of about 600 dbar.We concentrate on the North Atlantic instead of the global ocean because it is easier to see differences on a smaller scale and neutral physics are interesting in the North Atlantic due to the Mediterranean outflow producing increased values of neutral helicity.Showing  δz, will not be the final solution.We thus linearly interpolate S and onto the new surface and calculate new lateral density gradient errors.We then treat these new density gradient errors as we did the init -field before to get a more accurate optimized approximately neutral surface.If we repeat these steps often enough | | 2 will converge.Once it has converged the surface will be as close to neutral as possible, with the residual fictitious diapycnal diffusivity being due only to the path-dependency caused by neutral helicity.We will call this surface the ω-surface.
On all density surfaces we will have regions where the surface outcrops or hits the bottom topography.Due to this we will end up with several regions on a density surface which do not communicate with each other.A typical example is a marginal sea with narrow connections to the open ocean.In the algorithm described above we deal with this problem of independent regions by writing a set of equations as in Eq. ( 14) for each seperate region.Similar to before we constrain the average pertubation density of each region, , to be zero.
Ocean Sci., 5, 155-172, 2009 www.ocean-sci.net/5/155/2009/When optimizing approximately neutral surfaces with the method above we sometimes get a result where | | 2 does not converge.This is because of the algorithm overestimating the depth change, δz, due to the algorithm not knowing about the stratification above and below the surface optimized.This can then lead to a growing | | 2 due to the algorithm trying to overcorrect at these casts.If this happens we have to dampen the depth change; this means we only use a certain percentage of the depth change estimated by the algorithm to calculate the optimized surface.Another way of minimizing the possibility of this problem is to discard the data in the mixed layer -a region in which other processes than neutral physics are dominant.On all the following surfaces we will discard data shallower than 200 dbar.
A similar approach as above can be used to minimize for the slope error, s, instead of the density gradient error, .This would be more consistent with the aim of minimizing the fictitious diapycnal diffusivity but on the other hand the minimisation of the density gradient error is easier to understand when compared to the theoretical ideas in Appendix D. Both approaches give very similar resutls.
McDougall and Jackett (1988) contains an algorithm that similarly modifies existing approximate neutral surfaces by minimizing the size of the square of the density gradient errors, , at each spatial location, in this case weighted by N −2 .This was achieved using a multi-dimensional Newton technique, one dimension for each data point on the approximate surface, with one additional dimension for a Lagrangianmultiplier equation constraining the mean pressure perturbation to be zero.The computational method described above is a two-dimensional analogue of a new sparse matrix inversion technique that labels three-dimensional oceanographic data with a new neutral density variable γ i (Jackett et al., 2009).The optimization methods described in this paper and in McDougall and Jackett (1988) and Jackett et al. (2009) all have as their goal the minimization of (weighted) sums of squares • , the differences between the three methods being in the simplicity of the equations that are actually used.Mc-Dougall and Jackett (1988) used the set of linear equations to minimize • while assuming given values of the vertical gradients of salinity and potential temperature.The solution technique proceeded iteratively until convergence with revised values of the vertical gradients of salinity and potential temperature being made after each iteration if required.By contrast, the method of the present paper finds values of a logarithmic density perturbation, , such that the resulting • is minimized on the original surface in space.We then use this perturbation logarithmic density to estimate the pressure perturbation, as described by Eq. ( 16) above.This new surface is then iterated through the same process again until convergence is achieved.This description shows that the two methods are quite similar.We have found the present method to have good convergence properties and the code has been extended to include stations where the surface in question is not simply connected.As will be shown later, the de- Even though the change between our initial condition and the ω-surface are quite large in terms of the fictitious diapycnal diffusivity, the correlation between gN −2 H n and −∇ a ×ǫ•k or the variations of γ n on the ω-surface or the actual changes of temperature and pressure between the initial condition and the ω-surface are reasonably small.

How 'neutral' are existing density variables?
To show the differences between different density variables we use a density surface in the North Atlantic with an average depth of about 600 dbar.We concentrate on the North Atlantic instead of the global ocean because it is easier to see differences on a smaller scale and neutral physics are interesting in the North Atlantic due to the Mediterranean outflow producing increased values of neutral helicity.Showing  velopment of the optimization technique for a single surface leads to significant improvements in the accuracies achieved by the two-dimensional surfaces when compared with isosurfaces of three-dimensional variables (e.g. the code developed by Jackett and McDougall (1997) to calculate γ n ), all www.ocean-sci.net/5/155/2009/Ocean Sci., 5, 155-172, 2009 in terms of their abilities in approximating neutral tangent planes.The technique described here can be seen as a Lagrangian method, calculating the change of pressure of a surface, whereas in techniques used to label a three-dimensional data set can be seen as an Eulerian method where density is calculated at every point in x/y/z-space.
To illustrate the improvements of the optimized approximately neutral surface, the ω-surface, we choose a surface with an average pressure of about 1400 dbar.Pressure and conservative temperature on this surface are shown in Fig. 4a  and b.The surface chosen here is just an arbitrary example of a density surface covering the global ocean and the results are very similar for surfaces that are denser or lighter than the surface shown.Neutral helicity on the same surface is shown in Fig. 4c.The regions of elevated values of neutral helicity are mainly concentrated in the Southern Ocean (especially in the regions of high eddy activity) and in the North Atlantic (close to where the surface outcrops and close to the Mediterranean outflow).This is where we would expect high values of neutral helicity due to strong gradients of pressure and temperature.
Comparing Fig. 5a and b one can see the improvement achieved by using the algorithm introduced in this paper compared with the γ n -surface (which was used as initial condition).Shown is the fictitious diapycnal diffusivity, D f , plotted as log 10 D f where the colour scale was chosen to make the comparison of both surfaces possible.Both the North Atlantic and the Southern Ocean have regions with a fictitious diapycnal diffusivity larger than 10 −5 m 2 s −1 on the γ n -surface and therefore exceeding the values measured in most regions in the ocean.These are the regions where most density variables produce large errors with the other regions of the global ocean usually being less problematic.Most other regions have fictitious diapycnal diffusivities of approximately 10 −7 m 2 s −1 .This has been reduced by a few orders of magnitude in the ω-surface, pushing all the fictitious diapycnal diffusivities significantly below the values measured in the ocean with the remaining errors located close to the outcropping regions.On the ω-surface there are no fictitious diapycnal diffusivities larger than 10 −5 m 2 s −1 with most regions having values smaller than 10 −10 m 2 s −1 which is insignificant compared to the values measured in the ocean.The higher slope errors close to the outcropping regions are caused by high values of ∇ a p (and ∇ a ) causing high values of neutral helicity (compare Figs. 4c and 5b).
Similar results can be seen in Fig. 6 for data calculated from the WOCE climatology (Gouretski and Koltermann, 2004), comparing D f on a γ n -surface and an ω-surface with an average pressure of about 1400 dbar.As in the model output the new algorithm leads to an improvement (a reduction) in D f , even though it is not as large as in the example using model output due to the code used to calculate γ n (Jackett and McDougall, 1997) being dependent on a reference data set which is based on a Levitus climatology (Levitus, 1982), which is closer to the climatology used here than the model   output.Another difference between the model output and the climatology is the patchiness of the fictitious diapycnal diffusivity in the climatology.This is due to the averaging of observational data done to construct climatologies.
The improvement of the γ n -surface can also be seen by looking at the fictitious diapycnal diffusivity, D f , on a frequency plot (see Fig. 7).Fig. 7b shows the surface used in the previous figure and Fig. 7a and c show a lighter and a denser surface.One can see that D f decreases by a few orders of magnitude when using the ω-surface compared to the γ n -surface.The large improvement is possible because the model output we are using had water masses that deviated significantly from observed ocean properties.When applied to atlas data the fictitious diapycnal diffusivity in an ω-surface is perhaps just one to two orders of magnitude less than in a γ n -surface.
The improvement made by the algorithm can be seen by plotting γ n on an ω-surface (Fig. 8).γ n values have a range from 27.61 to 27.65 on the ω-surface which is a substantial density change.
Another way of seeing the improvement is by plotting T b ∇ a p×∇ a •k≈gN −2 H n vs. −∇ a × •k for the γ n and the ω-surfaces (see Fig. 9a and b for γ n and ω, respectively).For the ω-surface one can see a very good agreement between gN −2 H n and −∇ a × •k, as all the points of the surface almost end up on the line.
Even though the change between our initial condition and the ω-surface are quite large in terms of the fictitious diapycnal diffusivity, the correlation between gN −2 H n and −∇ a × •k or the variations of γ n on the ω-surface or the actual changes of temperature and pressure between the initial condition and the ω-surface are reasonably small.Ocean Sci., 5, 155-172, 2009 www.ocean-sci.net/5/155/2009/ 4 How "neutral" are existing density variables?
To show the differences between different density variables we use a density surface in the North Atlantic with an average depth of about 600 dbar.We concentrate on the North Atlantic instead of the global ocean because it is easier to see differences on a smaller scale and neutral physics are interesting in the North Atlantic due to the Mediterranean outflow producing increased values of neutral helicity.Showing only the North Atlantic also gives us the opportunity to use the density variable γ EW of Eden and Willebrand (1999), a density variable fitted only to the North Atlantic.The results shown below are very similar for depth ranges different to the average pressure of the density surfaces of about 600 dbar.
The density surfaces which we compare are the new and the old neutral density variables (γ i (Jackett et al., 2009) and γ n (Jackett and McDougall, 1997), respectively), a γvariable approximated with a rational function of salinity and conservative temperature (γ rf , McDougall and Jackett, 2005b), a γ -variable approximated with a function fitted to data of the North Atlantic (γ EW , Eden and Willebrand, 1999), potential density with reference pressures of 0, 600, 1000 and 2000 dbar (σ 0 , σ 600 , σ 1000 and σ 2000 ), orthobaric density (ρ ν , de Szoeke et al., 2000), modified steric anomaly surfaces and the optimized approximately neutral density surface, the ω-surface, of this paper.Note that the algorithm producing ω-surfaces improves exisiting density surfaces (i.e. it works in two dimensions), whereas all other density variables mentioned above label a three-dimensional hydrography with density.

Different approximations to neutral surfaces
Five different approximations to neutral surfaces have been discussed to date.All of them except γ EW are constructed to minimize s 2 or 2 (i.e.minimize the slope difference or density gradient errors between the approximately neutral density surface and the neutral tangent plane).
The first γ -variable, γ n , is dependent on a pre-labelled dataset and therefore the quality of a surface calculated with this technique is highly dependent on the proximity of the dataset to be labelled to the reference dataset (which is the Levitus climatology (Levitus, 1982)).Therefore if this code is used for model output simulating a different ocean (a paleo ocean or future climate) or if the model drifts from its initial state, the γ n variable may be less neutral than a well chosen potential density surface.This problem has been adressed with a new method of constructing approximate neutral surfaces, γ i , which uses the old γ n variable as an initial condition and an iterative inversion method to improve the surfaces.This new variable is computationally more expensive but significantly improves the accuracy of the surfaces.The third γ -variable, γ rf , is a rational function approximating neutral density surfaces dependent only on S and .In contrast to γ i and γ n , γ rf is independent of pressure, lat-   itude and longitude.Not being dependent on latitude and longitude means that it ignores the hemispheric changes in water-mass characteristics, therefore making it less neutral than the other neutral density variables (at least when used for a global density surface).The advantage of γ rf is that it is faster and easier to compute making it better for use by the ocean modelling community.γ EW is a neutral density variable constructed for use in the North Atlantic.Compared to the other approximate surfaces its main aim is not only to have the approximately neutral surface as neutral as possible but also to approximately satisfy the points mentioned in Sect.3, i.e. trying to make the horizontal gradient of the neutral density agree with the gradient of the in situ density and trying to make the vertical gradient of the neutral density proportional to the static stablility of the water column.Compared to the other γ variables, ω (as described in this paper) only improves a single surface rather than producing a continuum of surfaces in a three-dimensional dataset.
The γ i -surface (Fig. 10a) , which is the most accurate  method to date of achieving the neutral property in threedimensional hydrography, shows the smallest values of D f compared to all the other density variables in this analysis (apart from individual ω-surfaces constructed by the algorithm described in this paper).The main regions of increased fictitious diapycnal diffusivity are the Mediterranean outflow, the outcropping regions in the north and the Gulf stream region.We would expect increased values of slope error on a good approximate neutral surface in regions where we have increased values of neutral helicity, which is proportional to ∇ a p×∇ a .Such large values of neutral helicity are likely to occur in regions of either a strong pressure gradient on the surface, a strong temperature gradient or both.To have zero neutral helicity on a surface the pressure gradient and the temperature gradient have to be exactly aligned, or one of the two needs to be zero.The highest values of D f occur near Spain where there are strong pressure gradients and temperature gradients that do not align.Further off the coast of Spain there is a strong temperature gradient but the pressure gradient is quite small (the density surface is relatively flat) and therefore D f reduces drastically.The other two regions of high D f are mainly due to a very strong pressure gradient, near the outcropping of the density surface.
Looking at the γ n -surface (Fig. 10b) we can see a very similar pattern to the γ i -surface and slightly increased values of fictitious diapycnal diffusivity.These increased values are due to the offset of the model output from the reference data set used by the γ n -code as explained by McDougall and Jackett (2005b).This is the major problem of this density variable which has been adressed with the new γ i -variable Ocean Sci., 5, 155-172, 2009 www.ocean-sci.net/5/155/2009/The surface chosen for these plots has an average pressure of approx.600 dbar.The colours chosen for these plots are the same as for those of Fig. 10.(Jackett et al., 2009).The main improvements of γ i compared with γ n are in the Southern Ocean (not shown here) and the North Atlantic.
γ rf (Fig. 10c) gives a very small fictitious diapycnal diffusivity over most of the North Atlantic with the larger errors located at a concentrated region where the surface outcrops.This is likely due to a change in the outcropping region from the hydrography which has been used to construct this variable.
γ EW (Fig. 10d) is the γ -variable with the largest fictitious diapycnal diffusivity.The order of magnitude of this diffusivity is comparable with that of a potential density surface with a reference pressure which is not well chosen.The reason for this is that instead of trying to minimize only s 2 as with the other γ -variables, the aim of this function was to also minimize the other two mutually inconsistent points mentioned in Sect.3. ω (Fig. 10e) shows the smallest fictitious diapycnal diffusivity, with the highest values close to Spain.This twodimensional approach decreases this diffusivity by about two orders of magnitude, pushing D f far below the values measured in the ocean.The errors close to Spain are likely due to the crossland mixing scheme used in MOM4 to distribute the Mediterranean outflow into the North Atlantic.

Potential density
Potential density is a widely used density variable.At its reference pressure a potential density surface coincides with the neutral tangent plane but as soon as a potential density surface departs from its reference pressure, the slope of this surface increasingly differs from the slope of the neutral tangent plane.This can be seen by looking at the normal to the potential density surface, and the normal to a neutral tangent plane, realising that these two expressions are equal only at the reference pressure p r .
It can also be shown that the variations of potential density (referenced to p r ) along a neutral tangent plane are given by (Jackett and McDougall, 1997;McDougall and Jackett, 2005b) www.ocean-sci.net/5/155/2009/Ocean Sci., 5, 155-172, 2009 The fictitious diapycnal diffusivity on potential density surfaces referenced to 0, 600, 1000 and 2000 dbar are shown in Fig. 11.These potential density surfaces show a large fictitious diapycnal diffusivity in the east where the warm water of the Mediterranean enters and at the northern outcrop where cold surface waters are reached.
The fictitious diapycnal diffusivity in these figures is also due to the offset of the pressure on these surfaces from the reference pressure (Eq.19).The σ 600 -surface is the closest we can get to the approximate neutral surface due to the reference pressure being optimally chosen, with larger errors in the other potential density surfaces (proportional to the distance of the reference pressure to the average pressure of that surface).

Modified steric anomaly surfaces
A similar variable to potential density that has not been used much recently is steric anomaly (also called specific volume anomaly).Here we define a modified steric anomaly variable as where S r and r are fixed values of salinity and conservative temperature.This differs from the normal definition of steric anomaly by simply replacing 35 psu with some other fixed salinity and replacing in situ temperature of 0°C with a different reference temperature.It is important to note that once the reference parcel is decided on, the second part of the equation is a function only of pressure.
It can be shown that the variation of modified steric anomaly along a neutral tangent plane is given by where κ is the adiabatic and isohaline compressibility of seawater and ρ r and κ r are the density and compressibility at (S r , r , p) (where S r and r have been optimally chosen to minimize the spatial variation of δ).
The fictitious diapycnal diffusivity on the modified steric anomaly surface (Fig. 12a) shows the largest values close to large pressure gradients on the surface, which is a logical consquence of Eq. ( 21).The largest error is west of Spain close to where the crossland mixing scheme of MOM4 distributes the Mediterranean outflow into the North Atlantic.This region of large error can be seen in the fictitious diapycnal diffusivity of most surfaces but it is largest in the modified steric anomaly surface.The other region of large D f is along the highest ∇ a p -the region where the surface outcrops in the northern North Atlantic.One big advantage of using steric anomaly surfaces is the existence of a geostrophic streamfunction, the Montgomery potential (Montgomery, 1937).

Orthobaric density
Orthobaric density (de Szoeke et al., 2000), ρ ν , has recently been introduced as a density variable that is a function of pressure and in situ density that has the property that as long as water mass variations occur in a monotonic way with pressure along the neutral directions, it can be made quite neutral for a single ocean basin (McDougall and Jackett, 2005a).If used for the global ocean it is not possible to tune this variable so that it is a good approximation to neutrality.This is due to the inability of the variable to accurately accommodate differences between water masses at fixed values of pressure and in situ density such as occur between the Northern and Southern Hemisphere portions of the World Ocean (McDougall and Jackett, 2005a).
The fictitious diapycnal diffusivity for ρ ν is shown in Fig. 12c.Most of the large values for D f are concentrated at regions of highest ∇ ρ ν p.This can be seen by looking at the change of orthobaric density along a neutral tangent plane (Eq.14 of McDougall and Jackett, 2005a): where 0 is a reference conservative temperature and is an integrating factor, both a function of pressure and in situ density .ρ ν shows some of the largest errors in D f of all the surfaces analysed.

Further comparisons
Above we have seen two-dimensional maps of the fictitious diapycnal diffusivities on the previously described approximate density surfaces in the North Atlantic, giving us a view as to how good these surfaces are in representing isopycnal flow in the ocean.To further facilitate this intercomparison we now look at frequency distributions of D f .These distributions for γ n , σ 0 and ω, plotted on a log 10 -scale, can be seen in Fig. 13 .Fig. 14 shows the 95th-percentiles of the fictitious diapycnal diffusivities for all density surfaces previously considered.That is, the vertical axis shows the value of the fictitious diapycnal diffusivity of density that is exceeded by 5% of the data.
It is clear from Fig. 13 that there is a substantial decrease of fictitious diapycnal diffusivity going from σ 0 to γ n to ω.The surfaces on which these frequency distributions are calculated are the same as the density surfaces shown previously, all with an average pressure of 600 dbar.It is known that σ 0 is not the potential density surface with an ideally chosen reference pressure; by choosing a better reference pressure the fictitious diapycnal diffusivity would decrease, but still be larger than for the other surfaces shown.The γ nsurface is dependent on a reference dataset which limits its neutrality since the model output drifted from the reference data set.
Ocean Sci., 5, 155-172, 2009 www.ocean-sci.net/5/155/2009/where Θ 0 is a reference conservative temperature and Φ is an integrating factor, both a function of pressure and in-situ density .ρ ν shows some of the largest errors in D f of all the surfaces analysed.

Further comparisons
Above we have seen two-dimensional maps of the fictitious diapycnal diffusivities on the previously described approximate density surfaces in the North Atlantic, giving us a view as to how good these surfaces are in representing isopycnal flow in the ocean.To further facilitate this intercomparison we now look at frequency distributions of D f .These distributions for γ n , σ 0 and ω, plotted on a log 10 -scale, can be seen in Fig. 13 .Fig. 14 shows the 95 th -percentiles of the fictitious diapycnal diffusivities for all density surfaces previously considered.That is, the vertical axis shows the value of the fictitious diapycnal diffusivity of density that is exceeded by 5% of the data.The black vertical line shows a ficititious diapycnal diffusivity of 10 −5 m 2 s −2 ; values right of this line are larger than the mean value for the diapycnal diffusivity measured in the ocean.These values for the fictitious diapycnal diffusivity are for one actual surface (with an average pressure of 600 dbar).
It is clear from Fig. 13 that there is a substantial decrease of fictitious diapycnal diffusivity going from σ 0 to γ n to ω.The surfaces on which these frequency distributions are calculated are the same as the density surfaces shown previously, all with an average pressure of 600 dbar.It is known that σ 0 is not the potential density surface with an ideally chosen reference pressure; by choosing a better reference pressure the fictitious diapycnal diffusivity would decrease, but still Fig. 13.log 10 (D f ) of γ n (green), σ 0 (red) and ω (black).The black vertical line shows a ficititious diapycnal diffusivity of 10 −5 m 2 s −2 ; values right of this line are larger than the mean for the diapycnal diffusivity measured in the ocean.These values for the fictitious diapycnal diffusivity are for one actual surface (with an average pressure of 600 dbar).
In Fig. 14 one can see that for the North Atlantic the σ 2000 and ρ ν -surfaces have fictitious diapycnal diffusivities exceeding 10 −5 m 2 s −1 over more than 5% of their area.The decrease in fictitious diapycnal diffusivities from ρ ν -surfaces to potential density surfaces to approximate neutral density surfaces shows the considerable improvement that can be achieved by using more accurate density variables for inverse models or other applications of density surfaces describing isopycnal flow.

Conclusions
We have developed a new method for finding an individual approximately neutral surface through a three-dimensional hydrographic data set (either observational data or model output).The degree of non-neutrality along such an ωsurface has been minimized and the fictitious diapycnal dif- These values for the fictitious diapycnal diffusivity are for one actual surface (with an average pressure of 600 dbar).
be larger than for the other surfaces shown.The γ n -surface is dependent on a reference dataset which limits its neutrality since the model output drifted from the reference data set.
In Fig. 14 one can see that for the North Atlantic the σ 2000 and ρ ν -surfaces have fictitious diapycnal diffusivities exceeding 10 −5 m 2 s −1 over more than 5% of their area.The decrease in fictitious diapycnal diffusivities from ρ ν -surfaces to potential density surfaces to approximate neutral density surfaces shows the considerable improvement that can be achieved by using more accurate density variables for inverse models or other applications of density surfaces describing isopycnal flow.

Conclusions
We have developed a new method for finding an individual approximately neutral surface through a three-dimensional hydrographic data set (either observational data or model output).The degree of non-neutrality along such an ωsurface has been minimized and the fictitious diapycnal diffusivity associated with these ω-surfaces is the least that has been found to date using other surfaces.

The sm trality is inherent
The a and the surfaces ware to f au/TEOS Fig. 14.The 95th-percentiles of the fictitious diapycnal diffusivity in the North Atlantic is shown for all the density surfaces considered.One can see that for the North Atlantic the σ 2000 and ρ ν -surfaces have fictitious diapycnal diffusivities exceeding 10 −5 m 2 s −1 over more than 5% of their area.These values for the fictitious diapycnal diffusivity are for one actual surface (with an average pressure of 600 dbar).fusivity associated with these ω-surfaces is the least that has been found to date using other surfaces.
These surfaces are ideal to use as water mass density boundaries in inverse models.The improved neutrality of these surfaces might not be significant for large box inverse models of non-synoptic hydrographic sections but will definitely improve inverse models using synoptic sections for process studies that particularly target the determination of mixing.
The small deviation of these ω-surfaces from exact neutrality is shown to be limited by the neutral helicity that is inherent in the hydrographic data.
The algorithm for forming these ω-surfcaes is described and the extent of the non-neutrality of many other density surfaces is compared with these ω-surfaces.MAT-LAB software to form these surfaces is available at http://www.TEOS-10.org.www.ocean-sci.net/5/155/2009/Ocean Sci., 5, 155-172, 2009 where the line integral is in a surface A and the area element is normal to A, dA=m a dxdy.So if we can show the result will be proven (note that m a =k−∇ a z and is two dimensional so that •dl= •dr).
From McDougall and Jackett (1988) Eqs. ( 4) and ( 5) we have So the left-hand side of Eq. ( C3) is The first part of the right-hand side simplifies to ∇ a × •k because ∇ a × is an exactly vertical vector.Therefore our result C1 is proven.

Appendix D Equation 11 in detail
Here we derive a close connection between neutral helicity and ∇ a × .
Using the definition of (Eq.10) we get The first two terms on the right-hand side add up to zero because α S =−β and therefore ∇ a × becomes (using (D2) If we now use ∇ a S= α β ∇ a + 1 β from Eq. ( 10) we get The second term is small and will therefore be ignored from here on.( (from Eq. 9) so that Eq. (D3) establishes the desired approximate relation −∇ a × ≈gN −2 H n .Now we continue our derivation of Eq. ( 11) to result in T b A pd : From the first to the second line we use Stokes' theorem.The only approximation in Eq. ( D4) is the assumption that T b is constant so it can be taken outside the integral.

Appendix E Theoretical thoughts about the algorithm and numerical testcases
Here we explore the theoretical ideas on the relationship between neutral helicity and the density gradient errors which led to the development of the algorithm used to optimize approximately neutral surfaces.
We know that the path-dependent uncertainty does not occur on a single ocean section because we can link up all the neutral tangent planes on a section without any slope errors.Therefore we take a N-S section of the ocean and repeat it to the east and west so that there are no zonal gradients.This gives us an initial three-dimensional data set from which "perfect" neutral surfaces with neutral helicity being zero everywhere can be found.Having found a neutral surface from this artificial data set we can perturb a single point on this surface in a way that the perturbed bottle talks neutrally to the original water parcel, but neutral helicity is introduced.
An even simpler way to think of this numerical test case is to imagine a region of the ocean where the pressure and conservative temperature gradients on a particular approximately neutral surface are aligned so that ∇ a p and ∇ a are parallel everywhere.Now we consider a circular pertubation of in the presence of the backgound pressure gradient as shown in Fig. E1a.This pertubation--field produces a dipole of neutral helicity as shown in Fig. E1b.An isolated anomaly of or S, as shown in   If we now minimize 2 , we come to the solution seen in Fig. E3a (right), colour being the perturbed -field, , and the vectors showing new , the new field of density gradient errors with minimized 2 .Fig. E3a (left) shows ∇ a × .We find the curl at the right side of our block of density gradient errors has opposite sign to the curl at the left side.According to Eq. ( 11) we know that ∇ a × new is also approximately proportional to neutral helicity.This then resembles the neutral helicity dipole as seen in Fig. E1b, with both poles of this dipole moved apart more and with density gradient errors going around the neutral helicity poles.
To further look at the way the least-square inversion changes the density gradient errors we construct another test case.Fig. E3b shows a block of density gradient errors in which the left half of the box has initial density gradient errors of strength 12 and the right half has the same strength but opposite sign.As in the previous example we can see that ∇ a × new has high values along the right and left side of our blocks of imposed density errors, which corresponds to neutral helicity in the real ocean.The main change to the previous example is the border between the block of density errors with strength 12 and the block of density errors of strength −12, which gives us a ∇ a × new twice the strength of the other borders.The final density gradient vectors mainly circulate inside this region between the maximum/minimum values of ∇ a × new .The new vectors are spread beyond the initial white box and rotate around the strips of positive and negative helicity.

Fig. 1 :
Fig.1: This figure shows the differences between e, e hel , e a and w.The lateral velocity (ui + vj) is directed horizontally.The three surfaces shown are the approximately neutral surface (γ a ), the neutral tangent plane (ntp) and the top-most is the lateral velocity plus w, where w includes all components leading to a flow which differs from a purely horizontal flow (the tilt of an approximately neutral surface, mixing effects and a diapycnal velocity caused by the ill-defined nature of neutral surfaces, e hel ).

Fig. 2 :Fig. 3 :
Fig.2: The blue surface shows the initial approximately neutral surface (γ a -surface) on which the density pertubation field, Φ ′ , is calculated.This density pertubation field is converted into a depth change, δz, which is then applied to the initial surface to get the new apaproximately neutral surface (the blue surface).
The grid used by the algorithm explained in this paer.The red points are the tracer grid points and the slopes rrors/pressure gradient errors are calculated on the green oints.

Fig. 3 .
Fig.3.The grid used by the algorithm explained in this paper.The red points are the tracer grid points and the slopes errors/pressure gradient errors are calculated on the green points.
density for which | | 2 is minimized.To solve this set of equations, minimize|A − init | 2 , Fig. 5: log 10 (D f ) on the (a) γ n -surface and the (b) ωsurface.
for the γ n and the ω-surfaces (see Figs.9 (a) and (b) for γ n and ω respectively).For the ω-surface one can see a very good agreement between gN −2 H n and −∇ a × ǫ • k, as all the points of the surface almost end up on the line.
Fig. 6: log 10 (D f ) on the (a) γ n -surface and the (b) ω-surface calculated using WOCE climatology data.
log 10 (D f
Fig.7: Frequency plot of log 10 (D f ) for the γ n -surface (red) and the ω-surface (black) with an average pressure of (a) 1000 dbar, (b)1400 dbar (the surface used throughout the text) and (c) 1800 dbar.The black vertical line shows a ficititious diapycnal diffusivity of 10 −5 m 2 s −2 ; values right of this line are larger than the mean value for the diapycnal diffusivity measured in the ocean.

Fig. 7 .
Fig. 7.Frequency plot of log 10 (D f ) for the γ n -surface (red) and the ω-surface (black) with an average pressure of (a) 1000 dbar, (b) 1400 dbar (the surface used throughout the text) and (c) 1800 dbar.The black vertical line shows a ficititious diapycnal diffusivity of 10 −5 m 2 s −2 ; values right of this line are larger than the mean value for the diapycnal diffusivity measured in the ocean.
Fig. 8: γ n on the ω-surface for model output from a MOM4 model run.

Fig. 8 .
Fig. 8. γ n on the ω-surface for model output from a MOM4 model run.

Fig. 8 :
Fig. 8: γ n on the ω-surface for model output from a MOM4 model run.

Fig. 9 .
Fig. 9. T b ∇ a p×∇ a •k≈gN −2 H n vs. −∇ a × •k for (a) a γ n and (b) an ω-surface.The rms error of the difference between theory and the plotted data decreases by a factor of 6.

Fig. 10 .
Fig. 10.log 10 (D f ) for (a) γ i , (b) γ n , (c) γ rf , (d) γ EW and (e) ω.The surface chosen for these plots has an average pressure of approx.600 dbar.The same colour scale is used in each plot.

Fig. 11 .
Fig. 11.log 10 (D f ) for (a) σ 0 , (b) σ 600 , (c) σ 1000 and (d) σ 2000 .The surface chosen for these plots has an average pressure of approx.600 dbar.The colours chosen for these plots are the same as for those of Fig. 10.

Fig. 12 .
Fig. 12. log 10 (D f ) for (a) δ and (b) ρ ν .The colours chosen for these plots are the same as for those of Fig. 10 and Fig. 11.

Fig. 13 :
Fig.13: log 10 (D f ) of γ n (green), σ 0 (red) and ω (black).The black vertical line shows a ficititious diapycnal diffusivity of 10 −5 m 2 s −2 ; values right of this line are larger than the mean value for the diapycnal diffusivity measured in the ocean.These values for the fictitious diapycnal diffusivity are for one actual surface (with an average pressure of 600 dbar).

AFig. 14 :
Fig.14:The 95 th -percentiles of the fictitious diapycnal diffusivity in the North Atlantic is shown for all the density surfaces considered.One can see that for the North Atlantic the σ 2000 and ρ ν -surfaces have fictitious diapycnal diffusivities exceeding 10 −5 m 2 s −1 over more than 5% of their area.These values for the fictitious diapycnal diffusivity are for one actual surface (with an average pressure of 600 dbar).
Fig. E1: (a) A circular pertubation of Θ is shown in the presence of a background pressure gradient (the pressure gradient not shown; it increases linearly along the y-axis).(b) The resulting dipole of neutral helicity.

Fig. E3 :
Fig.E3: Solutions for the simplified density gradient errors as described in the text.

Fig. E3 .
Fig. E3.Solutions for the simplified density gradient errors as described in the text.
a γ n and (b) an ω-surface.The rms error of the difference between theory and the plotted data decreases by a factor of 6.
p≈∇ n p and ∇ a ≈∇ n so that T b ∇ a p×∇ a •k≈T b ∇ n p×∇ n •k=gN −2 H n α and | | is much less than |α ∇ a |). a