Generalized Patched Potential Density and Thermodynamic Neutral Density: Two New Physically Based Quasi-Neutral Density Variables for Ocean Water Masses Analyses and Circulation Studies

In this paper, two new quasi-neutral density variables—generalized patched potential density (GPPD) and thermodynamic neutral density g T —are introduced, which are showed to approximate Jackett and McDougall empirical neutral density g n signiﬁcantly better than the quasi-material rational polynomial approximation g a previously introduced by McDougall and Jackett. In contrast to g n , g T is easily and efﬁciently computed for arbitrary climatologies oftemperature and salinity (bothrealisticand idealized), has a clear physical basis rooted in the theory of available potential energy, and does not suffer from nonmaterial effects that make g n so difﬁcult to use in water masses analysis. In addition, g T is also signiﬁcantly more neutral than all known quasi-material density variables, such as s 2 , while remaining less neutral than g n . Because unlike g n , g T is mathematically explicit, it can be used for theoretical as well as observational studies, as well as a generalized vertical coordinate in isopycnal models of the ocean circulation. On the downside, g T exhibits inversions and degraded neutrality in thepolarregions,wheretheLorenzreferencestateisthefurthestawayfromthe actualstate.Therefore,while g T represents progress over previous approaches, further work is still needed to determine whether its polar de-ﬁciencies can be corrected, an essential requirement for g T to be useful in Southern Ocean studies, for instance.


Introduction
The problem of how best to construct a quasi-neutral density variable suitably corrected for pressure is a longstanding fundamental issue in oceanography whose answer is vital for many key applications ranging from the study of mixing to ocean climate studies.These include but are not limited to the separation of mixing into ''isopycnal'' and ''diapycnal'' components necessary for the construction of rotated diffusion tensors in numerical ocean models (Redi 1982;Griffies 2004), the construction of climatological datasets for temperature and salinity devoid of spurious water masses (Lozier et al. 1994), the construction of inverse models of the ocean circulation (Wunsch 1996), the tracking and analysis of water masses (Montgomery 1938;Walin 1982), the construction of isopycnal models of the ocean based on generalized coordinate systems (Griffies et al. 2000;de Szoeke 2000), the study of the residual circulation (Wolfe 2014), and the parameterization of mesoscale eddy-induced mass fluxes (Gent et al. 1995).
Physically, it is generally agreed that a suitable density variable g should possess the desirable dual thermodynamic, and dynamical attributes of defining adiabatic surfaces (the thermodynamic attribute) along which fluid parcels experience no net buoyancy force (the dynamical attribute; e.g., McDougall 1987;de Szoeke and Springer 2000;Huang 2014).The first attribute, which is equivalent to material conservation (also often referred to as quasimaterial conservation, meaning here conserved whenever u and S are also conserved), poses no difficulty as it can always be enforced by requiring g to be a function of potential temperature u and salinity S only.The second attribute-usually referred to as the neutral property-is problematic, however, as it can only be satisfied in special circumstances not usually encountered in the ocean.To satisfy exact neutrality, =g would need to be parallel at every point to the local neutral vector d 5 g(a=u 2 b=S) 5 2(g/r)(=r 2 c 22 s =p), where a and b are the thermal expansion and haline contraction coefficients defined relative to u, S, and p; g is the acceleration of gravity; c 2 s is the squared speed of sound; r is in situ density; and p is pressure.To understand why the latter property cannot be satisfied in general, it is useful to decompose =g into components parallel and orthogonal to d as follows: where b is an integrating factor, and R is a residual term perpendicular to d. Taking the curl of (1) and multiplying the result by d gets rid of =g and yields an equation for the residual R, namely, where the term H 5 d Á (= 3 d) is the helicity of the neutral vector d, which shows that exact neutrality can only be achieved when H 5 0, a well-known result (McDougall 1987;de Szoeke and Springer 2000;Huang 2014), with Eden and Willebrand (1999) discussing some of the conditions necessary for the helicity to vanish.In practice, achieving H 5 0 in the ocean would either require the ocean to be at rest-as r would then be a function of pressure p alone-or in absence of densitycompensated temperature-salinity variations along surfaces g 5 constant, which is equivalent to say that the ocean would then have a well-defined temperaturesalinity relationship of the form u 5 u(g) and S 5 S(g).
In the ocean, however, the existence of densitycompensated u/S anomalies conspire with thermobaricity (the pressure dependence of the thermal expansion coefficient) to make H nonzero and hence forbid the construction of exactly neutral density variables.McDougall and Jackett (1988) discuss the way the ocean appears to conspire to keep values of helicity low.
If so, what then are the physical principles determining the degree of nonneutrality that g should have?In particular, should material conservation be retained, or sacrificed to improve neutrality?From a theoretical viewpoint, the natural starting point for constructing a density variable are the thermodynamic equations for salinity and potential temperature (or Absolute Salinity and Conservative Temperature) combined into the following equation for density: where r(S, u, p) is viewed as a function of salinity, potential temperature, and pressure; q represents changes in density due to molecular diffusive effects of heat and salt.The left-hand side of (3) defines the differential form d-5 dr 2 c 22 s dp, which in general is not perfect and hence not integrable because of the nonzero helicity of the neutral vector; otherwise, it is well accepted that -(possibly modified via the introduction of an appropriate integration factor) would define the most natural choice of quasi-neutral, pressure-corrected density variable.Mathematically, this is equivalent to state that the total differential dg of any mathematically well-defined quasi-neutral density variable g can at best be written in the form and involves a nonvanishing, residual, imperfect differential form dw, with b as an integrating factor.Equation (4) yields (1) upon making the following substitutions dg / =g, du / =u, dS / =S, and dw / R as well as by interpreting S and u as their climatological values rather than their instantaneous ones.The main advantage of (1) is that it defines the problem in standard Euclidean vector space, which makes it easy to define the ''smallness'' of the residual R or its orthogonality with the neutral vector d.Working with ( 4) is harder mathematically, as it is less easy to define a norm in the space of differential forms that can be used to say whether dw is small or large.
The aforementioned mathematical difficulty has so far prevented the discovery of the ''right'' way of integrating the density equation [( 4)], with standard potential density, patched potential density (PPD), and orthobaric density representing the most well-known attempts.None of these density variables, however, is regarded as fully satisfactory.In absence of any clear theoretical argument on how best to approach (4), McDougall's (1987) postulate that the best quasi-neutral density variable should be one constrained to be as neutral as feasible, which in practice can be constructed by means of the Jackett and McDougall (1997, hereinafter JMD97) neutral density software, has been widely accepted as the most appealing alternative. 1evertheless, g n has a number of shortcomings that tend to restrict its use primarily to observational studies of the present-day ocean outside such regions as the Arctic Ocean or Mediterranean Sea, where it is currently not defined.Indeed, its relatively high computational cost makes its use prohibitive in numerical ocean modeling studies; its lack of mathematically explicit form forbids its use in theoretical studies of the ocean circulation, and its nonmateriality makes its use in inverse studies of ocean mixing or in the analysis of water masses using Walin's (1982) approach conceptually problematic, owing to the difficulty of evaluating thermobaric dianeutral dispersion rigorously.Moreover, the physical basis for g n remains arguably quite unclear.Indeed, density variables-such as potential density or orthobaric density, for instancetend to be defined as purely thermodynamic concepts having (or not) desirable dynamical properties when used for recasting the equations of motion in thermodynamic coordinates, for example, de Szoeke (2000).Thus, de Szoeke and Springer's (2000) orthobaric density is defined as purely thermodynamic variable function of in situ density and pressure only, which dynamically defines an exact geostrophic streamfunction.In contrast, JMD97's construction of g n tends to emphasize (approximate) dynamics over thermodynamics by defining it primarily in terms of the equation b 5 2d Á dx ' 0, where b is meant to represent the buoyancy of a single parcel experiencing adiabatic and isohaline lateral displacements dx.The thermodynamic properties of g n remain unclear, however, and a controversial topic [e.g., McDougall andJackett (2005b) vs de Szoeke andSpringer (2009)].One of several difficulties is that neutral trajectories obtained as solutions of d Á dx 5 0 implicitly require nonmaterial sources of heat and salt, thus violating the assumptions underlying the definition of buoyancy b 5 2d Á dx, since, as showed by McDougall (1987), they generally end up at a different vertical position than they originate from if integrated over a closed loop (around the main gyre of the North Atlantic ocean for instance).
Although JMD97 chose not to enforce materiality as a way to maximize the neutrality of their variable g n , both McDougall and Jackett (2005b) and de Szoeke and Springer (2009) seem to agree that nonmateriality is an undesirable feature of a quasi-neutral density variable, since it may confound the determination of diapycnal mixing in inverse ocean modeling studies for instance.So far, however, while it is generally agreed that a purely material function of S and u can be constructed to be quite neutral over a limited region of the ocean, as showed by Eden and Willebrand (1999) for the North Atlantic Ocean, McDougall and Jackett (2005b) have speculated that this is fundamentally impossible to achieve for the global ocean, after failing to construct a quasi-material rational polynomial approximation g a of g n , whose neutrality appeared to be no better than that of s 2 , while being neither a good approximation of g n , nor of its gradient.
The main novelty of the present paper is to show that JMD97 empirical neutral density g n , despite being primarily based on heuristic considerations, actually contains useful information about how best to integrate (4).This is shown here by showing that g n is very close to a physically based quasi-neutral density variable that outperforms all known density variables in terms of neutrality.This variable is also materially conserved and naturally approximates g n significantly better than g a .Such a variable is called thermodynamic neutral density and is a function of the Lorenz neutral density that enters the theory of available potential energy, whose construction for a realistic ocean with a fully nonlinear equation of state was recently discussed by Saenz et al. (2015).To that end, our paper proceeds in two steps: The first step, detailed in section 2, provides a new look at the concept of patched potential density, which is the concept that historically prompted the construction of neutral density and orthobaric density.This section argues that the classical expression for patched potential density is not a useful one for lacking any information about the actual patching process whereby density surfaces in different depth ranges are joined up at discontinuity points.An improved PPD, called generalized patched potential density (GPPD), which is significantly less discontinuous than the original PPD and which explicitly accounts for the patching process, is constructed.The advantage of GPPD is to make immediately clear what its continuous limit should be.Thermodynamic neutral density is one particular example of continuously differentiable analog of GPPD, whose construction, comparison with other density variables, and neutral properties are discussed in section 3. Section 4 discusses the results and their implications.

A generalized patched potential density
explicitly accounting for the patching process

a. Statement of the problem
Let us assume, like JMD97, de Szoeke and Springer (2000), and many others before, that the concept of patched potential density holds the key for understanding how best to construct a quasi-neutral density variable corrected for pressure.To that end, let us first recall that the standard definition of PPD is potential density referenced to a piecewise constant pressure field, namely, where i 5 1, . . ., N 2 1, where N is the number of discrete depth levels entering the construction of PPD.To be specific, let us assume z 1 5 0 and z N is a depth typical of the abyssal ocean.Here, we use the suffix prior to draw the attention of the reader to the fact that such a definition of PPD is merely a front window for what PPD is really about, given that (5) does not actually tell us anything about the patching process whereby density surfaces in different depth ranges are joined up at points of discontinuity.Although the patching process itself can be described in mathematical terms, it would be complicated and tedious to do so.For the present purposes, we adopt a much simpler approach to the patching process, which we believe captures its essence while being mathematically clearer.Our approach stems from the realization that the patching process is fundamentally about dealing with the discontinuous character of PPD prior , in the sense that if that PPD prior were continuous, one would just trace water mass properties along surfaces of constant PPD prior .As a result, we redefine the patching process as a process fundamentally aimed at making PPD prior as continuous as feasible through the successive removal of discontinuities, ultimately leading to the following posterior form of PPD: where p ijk and s ijk are possibly three-dimensional piecewise constant fields.We call (6) GPPD and regard it as the postpatching process (posterior) or true form of PPD.As we show below, a careful choice of the piecewise constant p ijk and s ijk can make GPPD significantly less discontinuous than PPD prior , thus enabling GPPD to appear reasonably continuous when plotted if enough discrete elements are chosen.
b. Description of the patching process as the successive removal of discontinuities Our approach to the patching process aims to make PPD prior as continuous as feasible without altering significantly =PPD prior , which is what makes PPD prior quasi neutral.To that end, we initiate the patching process by first subtracting a piecewise constant density offset in each depth range as follows: where p r (z) and s offset ( p) are assumed to be C ' functions of depth and pressure, respectively.The jump in PPD (1)  across a discontinuity, denoted [PPD (1) (z k , z k11 )], becomes where S b and u b are salinity and potential temperature values along the surface of discontinuity.The natural choice to reduce the discontinuity is to choose but unless S b and u b obey a well-defined u/S relationship of the form S b 5 S b ( p) and u b 5 u b ( p), the approach will only succeed at removing the discontinuity locally, not globally.To cure the problem, it is necessary to make the density offset vary with horizontal position as well, suggesting the following second modification: with x i , i 5 1, . . ., N i , y j , j 5 1, . . ., N j , a series of discrete points in the horizontal directions aimed at capturing the spatial variations in the u/S relationship.The problem, however, is that making s offset vary with horizontal position must in turn introduce horizontal discontinuities in PPD (2) , which can only be corrected by making p r vary with horizontal position as well, leading to the following (and last) modification: which is consistent with the form of generalized patched potential density [(6)] given above.Obviously, this is as far as the patching process can go, as we have run out of options for further curing discontinuities.

c. Validation
The above description of the patching process is arguably only qualitative.In practice, a full implementation of the method would require writing down explicit equations for the horizontal and vertical density jumps as well as providing an explicit procedure for constraining the number of discrete elements and the values of the piecewise constant p ijk and s ijk in (6).This is not further pursued here, however, as our primary aim is to use the concept of GPPD as a stepping stone for clarifying the continuous limit of PPD and introducing the concept of thermodynamic neutral density discussed in the next section.Before we do that, however, we first seek to validate the concept of GPPD.Specifically, if our hypothesis that (6) represents the true or revealed form of PPD, it should be possible to construct the piecewise constant fields p ijk and s ijk to obtain an accurate approximation of JMD97 empirical neutral density g n ; indeed, as shown by the latter authors, g n is known to behave as PPD; if so, g n should also behave like GPPD.The following aims to show that this is indeed the case and that a good agreement can, in fact, be achieved by choosing p ijk and s ijk to vary with latitude and depth only.
To that end, using the g n field supplied as part of the Gouretski and Koltermann (2004) WOCE dataset, we computed the 2D fields p ijk and s ijk , minimizing the misfit between g n and g GPPD using all possible data points for which g n is defined for an a priori-given partition V jk .Through trial and error, we settled on the particular twodimensional partition of the ocean volume depicted in the top-left panel of Fig. 2 (shown below), with Dz 5 500 m and Dy ' 208, using a least squares approach to find the optimal values of p jk and s jk in each subdomain.The main intent here is only to illustrate the feasibility of GPPD, not to explore systematically the sensitivity of the results to the different possible choices of volume partition or constraints on p ijk and s ijk .
The results are illustrated in Fig. 2. Interestingly, the top-left and bottom-left panels strongly suggest that s jk is primarily controlled by p jk at leading order, which is confirmed by the regression analysis depicted in the bottom-right panel.Although it would be in principle possible to modify our optimization procedure to constrain the piecewise pressure values p ijk to be close to actual pressures, this was not done here, in order to see whether the procedure would do it on its own or not.Rather, p ijk was just imposed to lie within the range of pressures encountered in the ocean.It is therefore interesting to see that rather than choosing a piecewise pressure field close to the reference pressure field p r (z k ) used in PPD prior and a density offset a function of both horizontal position and p r , the optimization procedure naturally chooses a pressure field that can depart occasionally strongly from p r (z k ), with a density offset function that is simply a one-dimensional function of the latter.The associated plot for g GPPD is given in the topleft panel of Fig. 1 for the 308W latitude-depth section in the Atlantic Ocean, which can be compared with the corresponding section for g n in the top-right panel.The strong similarity between the two figures is striking, given that the ability of g GPPD to reproduce the main features of g n is achieved with only 7 3 11 5 77 discrete reference pressures p jk ; the visual agreement is further confirmed by the scatterplot of g n against g GPPD depicted in the bottom panel of Fig. 4 (see below), which shows a near-perfect correlation between the two quantities [the outliers seemingly originating from somewhat strange values of WOCE g n in enclosed seas, for which the use of JMD97 neutral density software is a priori not valid].A histogram of the differences g GPPD 2 g n (blue bars in top panel of Fig. 4) shows that g GPPD approximates g n to better than 0.01 kg m 23 in most of the ocean, which is remarkable, as this is much better than what is achieved by McDougall and Jackett's (2005b) variable g a (red bars in top panel), which was specifically constructed to best approximate g n .
Intriguingly, the structure of the p jk 's is in fact much more reminiscent of that of the reference pressures that fluid parcels would have in their reference state of minimum potential energy that have been recently described in some recent advances in APE theory by Tailleux (2013b) and Saenz et al. (2015).The possibility to use APE theory to provide a physical basis for p r is confirmed in the next section and suggests that part of the structure of neutral density can be explained in terms of Lorenz reference density.

Continuous limit of PPD and connection with
Lorenz theory of available potential energy a. Implications for the continuous limit of PPD A classical result of analysis is that all continuous functions can be viewed as the limit of piecewise constant functions.It follows that by making the number of discrete elements forming the volume partition V ijk arbitrarily large, the piecewise constant pressure and density offsets can be assumed to converge toward continuous fields p ijk / p r (x) and s ijk / s r (x), respectively.This in turn implies GPPD / GPPD continuous 5 s[S, u, p r (x)] 2 s r (x) , ( 12) which is hence also continuous.In the following, we further constrain the form of the density offset by assuming it to be a function of p r alone, that is, s r (x) 5 s r1d [ p r (x)], as in the previous section.Mathematically, the continuous limit of =(GPPD) is not so well defined, since in general all what can be ascertained for the limit of a piecewise constant function is that it be continuous, not continuously differentiable (differentiability is associated with piecewise linear functions, not just piecewise constant).In other words, it is usually not true that =(GPPD) converges uniformly toward =(GPPD continuous ), where As shown below, this mathematical difficulty is reflected in the fact that the function g T introduced below does a better job at approximating g n than its gradient.
Motivated by the results of the previous section drawing a link between the structure of p ijk and that of the reference pressure entering Lorenz APE theory, we introduce a new quasi-neutral density variable, called thermodynamic neutral density g T , defined as where p LZ r (S, u) is the reference pressure that a parcel would have if brought in a notional reference state of rest obtained by means of an adiabatic and isohaline rearrangement of the actual state.As shown recently by Tailleux (2013b) and Saenz et al. (2015), the reference pressure p LZ r that fluid parcels would have in the Lorenz reference state of minimum potential energy r LZ r (z) is the solution of the level of neutral buoyancy (LNB) equation where the possible time dependence of the reference state, for example, Tailleux (2013a), is neglected for the moment.Future work, however, should aim to establish when and where the temporal variations of Lorenz reference state matter.For a time-independent reference state, the LNB equation [( 15)] implies that the reference depth of fluid parcels z r 5 z r (S, u) is a materially conserved quantity; solving (15) at all points in the ocean provides the following explicit construction for the continuous reference pressure field p r (x), namely, The reference density profile r LZ r (z) was estimated for the WOCE dataset following the methodology detailed in Saenz et al. (2015), with an example of the resulting p r (x) field at 308W in the Atlantic Ocean being illustrated in the top-right panel of Fig. 2.

b. Comparison between g T and g n
To compare g T with g n , we calibrated the unknown density offset function s r1d (p r ) to minimize the differences between the two variables, in order to make g T traceable to g n as defined in Huber et al. (2015).Traceability is important in order to interpret the remaining differences between the two variables as due to differences in the physics rather than due to some of the arbitrary choices that enter the construction of density variables.This was done here by means of a joint pdf analysis of the respective distributions of r(S, u, p LZ r ) and g n , with s r1d (p r ) constructed so as to minimize the misfit between g T and g n .At leading order, s r1d (p r ) is found to behave primarily as a linear function of z r by working in the context of the Boussinesq approximation, using p r 5 2r 0 gz r , with r 0 5 1035.40715 kg m 23 and g 5 9.81 m s 22 , the value of r 0 being obtained by regressing the pressure and depth fields provided as part of the WOCE dataset.After some trial and error, we finally settled on the following fit for s r1d (p r ): with a 5 20.023 364 812 862 605 and b 5 0.004 527 584 358 902, where the leading-order linear behavior of s r1d ( p r ) was further corrected by the piecewise polynomial function P oly depicted in Fig. 3.This piecewise polynomial is made up of two distinct polynomials, whose slopes at the intersection point very discontinuously.Attempts at removing this discontinuity degraded the agreement between g n and g T in the deep ocean, suggesting that g n might be functionally different in the deep ocean relative to the rest of the ocean, perhaps as the result of difficulties with the method of characteristics underlying its construction in presence of bottom topographic features.The distribution for g T obtained from such a procedure is depicted in the bottom left of Fig. 1 for the same Atlantic Ocean section at 308W previously used.An atlas providing plots of g T every 108 of longitude, along with the corresponding plots for g n and McDougall and Jackett (2005b) variable g a is also made available as supplemental material.Clearly, g T does much better than g a at capturing the main features and details of g n .The main differences primarily occur in the upper region of the ACC where g T displays some inversions not seen in g n .At the same time, it seems important to point out that as explained in JMD97, the values of g n in the Southern Ocean are not obtained from the actual Levitus data but from modified ones, which are found necessary for correctly labeling neutral densities of the second kind in such a region and which may explain some of differences between g T and g n seen there.Apart from this issue, Fig. 4 (bottom panel) shows that g T and g n are otherwise extremely well correlated.The upper panel shows a histogram of the differences between the two variables, which reveal that g T is, in general, better than g GPPD at approximating g n , although it also reveals a few instances of rather large differences between g T and g n that do not exist for g GPPD .
Another way to compare g T and g n is directly in (u, S) space.Although g n is not materially conserved, it is nevertheless possible to write it as a sum of a materially conserved part g n material (S, u) plus some residual dg.For the present purposes, we estimated g n material (S, u) as the bin average of g n in (u, S) space, using DS 5 0.1 psu and Du 5 0.18C for the binning, which is equivalent to defining g n material as the materially conserved function of u and S that best approximates g n in a least squares sense.Figure 5, top-left, top-right, and bottom-left panels show g n material , g T , and their residual, respectively.For comparison, the residual between g n material and g a is also provided on the bottom-right panel, which shows significantly higher values.Remarkably, g T and g n material appear to exhibit the same functional dependence on S and u for most of the ocean water masses, suggesting that the nonmateriality of g n might be the primary cause for the observed differences between g T and g n , even though the residual g T 2 g n appears to have a rather complex structure.Since the estimation of the nonmateriality of g n has proven so far technically complex and controversial [see de Szoeke and Springer (2009) vs McDougall and Jackett (2005b)], the present results are interesting as they might point to a potentially much simpler way to quantify the nonmateriality of g n , which we plan on investigating in future studies.

FIG. 4. (top)
Histogram of the decimal logarithm of the absolute value of the differences between g n and g GPPD (blue), between g n and g T (green), and between g n and g a (brown).(bottom) Scatterplots of g GPPD , g T , and g a against WOCE g n .The straight line is the 1:1 line of equation g n 5 g proxy , where the latter is either one of g T , g GPPD or g a .
To conclude this section, it is important to point out that the structure of the differences between g T and g n is somewhat sensitive to the way-by no means uniquethat the function s r1d ( p r ) is constructed and hence that these differences should only be regarded as indicative rather than definitive.Indeed, it is important to note that the WOCE g n is a neutral density of the second kind (NDSK) rather than of the first kind and not necessarily as neutral as could be.Moreover, there might be alternative ways to construct s r1d ( p r ) that would result in an even better agreement between g n and g T .On the other hand, it is also important to recognize that rather than constructing s r1d ( p r ) to minimize the differences between g T and g n , one might prefer to define it based on physical arguments.The most natural approach would be in terms of a globally defined u/S relationship parameterized in terms of p r , that is, of the form S r ( p r ) and u r ( p r ), which would yield This approach, however, is beyond the scope of the present paper and will be discussed in a subsequent study, along with the importance of retaining or not the temporal variations of Lorenz reference state.Note that an approach such as (18) becomes necessary when constructing thermodynamic neutral density for water mass distributions differing from present-day ones, since calibrating s r1d ( p r ) to mimic g n works only when g n is available (recall that the neutral density software is only valid for present-day climatologies).
c. How neutral is g T relative to other density variables?
Since JMD97's construction of neutral density focuses on =g n and its closeness to neutrality rather than on g n itself, it is of interest to examine to what extent the good agreement between the values of g T and g n demonstrated in the previous section also extends to the gradients.That this should be so is not mathematically guaranteed because it is easy to find counterexamples of two continuously differentiable functions, f(x) and g(x) being approximately equal to each other without this being true of their derivatives.2Moreover, as mentioned before, the process of defining the continuous limit of =(GPPD) is not as well defined as it is for GPPD itself.As a result, scatterplots of jj=g T jj versus jj=g n jj or of angle(=g T , d) versus angle(=g n , d) are expected to exhibit significantly more scatter than plots of g T versus g n , which was indeed verified.These are not shown as they are deemed to not be very informative, even though they tend to suggest that g T performs better than g a at approximating =g n .
A more telling diagnostic is illustrated in Fig. 6, whose aim is to illustrate the relative degree of neutrality of g T relative to that of s 0 , s 2 , g a , and g n , which is based on computing the data-based frequency distribution of the sine of the angle made between the gradient of a given density variable and the neutral vector: As expected (from Fig. 6d), the WOCE g n is the more neutral of all density variables, but the main new finding here is that g T is a surprisingly close second, clearly outperforming s 0 , s 2 , and g a , while achieving a degree of neutrality comparable in many ways to that of g n .Importantly, g T thus appears as the first globally defined materially conserved density variable to pass McDougall and Jackett's (2005b) ''kiss of death'' s 2 neutrality test, which they used to fail the de Szoeke and Springer (2000) orthobaric density.Perhaps more astonishing, however, is that g T naturally outperforms g a , a materially conserved density variable that was specifically constructed by McDougall and Jackett (2005b) to best approximate g n and its gradient.The result is important because it demonstrates that McDougall and Jackett (2005a,b) significantly underestimated the ability of a materially conserved variable to approximate neutral density.The question that the present findings raises, and which future studies should aim to clarify, is whether g T is the optimal way to construct a materially conserved approximation of g n or whether an even better way exists?
d.A posteriori rationalization of the relevance of Lorenz reference state to the theory of quasi-neutral density variables The strong agreement found between g n and g T confirms our hypothesis that the empirical neutral density procedure designed by McDougall (1987) contains important information about the physics of quasi-neutral density variables, which the present results suggest point to Lorenz APE theory.Can this be rationalized a posteriori?To see this, let us imagine that we have been able to find a solution g(S, u) to the neutral density equation constrained to be materially conserved, and let us show that it must necessarily be a function of Lorenz reference density.To that end, let us consider the set of all possible surfaces g(S, u) 5 constant; each of these surfaces can be plotted individually as depicted in the left panel of Fig. 7 and will in general have a complicated shape in the actual state of the ocean.Because all the surfaces g(S, u) 5 constant are materially conserved, each of the parcels making up such surfaces will remain on such surfaces in any adiabatic and isohaline rearrangements of the actual state, including the notional state of rest entering the Lorenz theory of available potential energy.We know, however, that in a state of rest, the isosurfaces of any solution to the neutral density equation must coincide with constant geopotential surfaces, as otherwise they would not be neutral.This implies therefore that g(S, u) must be a function of Lorenz reference density and hence that Lorenz APE theory is the natural way to think about quasi-neutral density variables if constrained to be materially conserved, thus confirming that the differences between in g n and g T must represent a measure of the nonmaterial conservation of g n .
An important caveat, however, is that the above proof is probably only valid as far as the natural vertical ordering of fluid parcels remains the same in the Lorenz reference state and the actual state.Because of thermobaricity, this vertical ordering of fluid parcels may occasionally undergo significant modifications in regions where the actual state becomes too far away from Lorenz reference state, causing g T to develop inversions, as is the case in the polar regions.Our hypothesized link between g n and g T , therefore, is likely to exist only in those regions of the ocean where the vertical gradients of the two quantities have both the same sign and both correctly predict ocean stability as measured by N 2 .

Discussion
In this paper, we have reexamined several issues pertaining to the construction of quasi-neutral density variables, starting with a fresh look at how best to define the continuous limit of patched potential density, a point of contention between JMD97 and de Szoeke and Springer (2000).To that end, it was found essential to redefine patched potential density because even though the latter is commonly defined as potential density referenced to a piecewise constant reference pressure field s[S, u, p r (z i )], z i # z # z i11 , this formula obscures the whole machinery of the actual patching process whereby different potential density surfaces are joined up at points of discontinuity.As a result, we introduced the concept of generalized patched potential density (GPPD) as the true or posterior form of PPD by reinterpreting the patching process as being primarily about the successive removal of the discontinuities of the prior form of PPD.Mathematically, GPPD relies on a threedimensional discrete partition of the total ocean volume and the specification of the three-dimensional piecewise constant pressure fields p ijk and s ijk .GPPD is significantly less discontinuous than PPD, so that it will look reasonably continuous when plotted.As a result, one may in principle trace water mass properties simply along surfaces of constant g GPPD , as one would naturally do with any standard continuous density variables such as s 2 .While practitioners of PPD may have a hard time reconciling their approach to the patching process with ours, or to accept that g GPPD really represents the true form of PPD, our approach is vindicated by showing that it is possible to construct g GPPD so that it approximates the WOCE g n dataset to better than 10 22 kg m 23 over most of the ocean, which is significantly better than what can be achieved by McDougall and Jackett's (2005b) quasi-material rational polynomial approximation g a .If g n behaves both like PPD and FIG. 7. Schematics of the argument establishing that neutral density needs to be a function of Lorenz reference density when constrained to be materially conserved.
GPPD, GPPD must represent a valid alternative way to look at PPD.
If one accepts our proposal to regard GPPD as the true form of PPD, then the problem of how best to define the continuous limit of PPD becomes trivial, the latter being then necessarily given by g GPPD continuous 5 s[S, u, p r (x)] 2 s offset (x) or g GPPD continuous 5 s[S, u, p r (x)] 2 s offset [x, y, p r (x)].Again, this conclusion is vindicated by the possibility of constructing a very good (continuously differentiable) approximation of the WOCE g n dataset by using as p r the reference position of a fluid parcel in Lorenz reference state, thus motivating the introduction of a new materially conserved quasi-neutral density variable, called thermodynamic neutral density, defined as where the function s r1d (p r ) was defined as a piecewise polynomial function of p r calibrated to minimize the differences with WOCE g n .The main finding here is that g T approximates g n and its gradient considerably better than g a (S, u), which was specifically designed to best approximate g n .But the ability of g T to mimic g n does not stop there.In addition to approximating g n very well, g T possesses shares additional features that until now were thought to be the sole prerogative of g n , namely, as follows: d The variable g T is the first physically based globally defined density variable that significantly outperforms all other density variables, apart from g n , in terms of neutrality; in particular, it passes the McDougall and Jackett (2005b) s 2 neutrality kiss of death test.
d The ACC region poses similar problems to both g T and g n .Thus, the ACC region is associated with the possibility of multiple neutral paths for g n (JMD97), whereas it is associated with multiple levels of neutral buoyancy for g T (Saenz et al. 2015).Moreover, it is where g T displays inversions not seen in g n , whereas it is where the Levitus dataset used by JMD97 needs to be modified in order for their neutral density software to function correctly.
d Like g n , g T possesses interhemispheric differences in water mass properties, so that it is not affected by what McDougall and Jackett (2005a) consider to be a challenge for quasi-neutral density variables.
d Like g n , g T is affected by both cabbeling and thermobaricity, in contrast to standard potential density, as discussed in Iudicone et al. (2008).
d As is well known, g n is a function of thermodynamic variables u, S, and p as well as horizontal position, whereas if the time dependence of Lorenz reference state is retained, g T (S, u, t) is also dependent on time, although not on space.Moreover, because of thermobaricity and the existence of multiple levels of neutral buoyancy in some parts of u/S space, it is in principle possible for the Lorenz reference state to change with time purely as the result of adiabatic changes, which could potentially cause adiabatic vertical dispersion with no signature in microstructure measurements, as is believed to be the case for g n (but for physically quite different reasons).
The fact that the physically based g T is naturally capable of approximating g n significantly more accurately than McDougall and Jackett's (2005b) rational polynomial approximation g a , while also sharing most of g n 's key attributes, strongly suggests that Lorenz reference state should play a more important role in the theory of neutral density than previously realized, at least in those regions of the ocean where the vertical gradients of the two quantities both have the same sign and both correctly predict ocean stability.This view appears to be supported, at least partly, by the fact that whereas neutral density has so far represented the main basis for thinking about how to define isopycnal and diapycnal directions in the ocean, it is the theory of available potential energy that has formed the main basis for the rigorous study of diapycnal mixing in the stratified turbulent mixing community, following the pioneering work of Winters et al. (1995).Moreover, while mesoscale eddy parameterizations generally rely on isopycnal directions based on the local neutral tangent plane, Gent et al.'s (1995) view that mesoscale eddies should act as a sink of APE suggests that such parameterizations might be more naturally formulated based on g T .On the other hand, it is important to point out that Winters et al.'s (1995) APE framework has been so far validated only for a linear equation of state, for which the concept of density is unambiguously and uniquely defined; in contrast, the number of quasi-neutral pressure-corrected material density variables of the form g(S, u) in a thermobaric ocean in the presence of density-compensated u/S anomalies is potentially infinite.Moreover, it can also be argued that it is the probability density function (PDF) attribute of the Lorenz reference state that is really the property that matters for diagnosing mixing rigorously rather than its connection to available potential energy.Thus, one could argue in the oceanic case that diagnosing mixing rigorously could equally well be achieved by analyzing the temporal behavior of the PDF of s 0 , s 2 , or s 4 [although how to relate the effective diffusivity in PDF space, e.g., Winters et al. (1995), to observed values of diapycnal mixing in physical space is a priori not straightforward].Moreover, while the vertical gradient of Lorenz reference density is in

FIG. 1 .
FIG. 1. (top left) GPPD based on the GPPD reference pressure and density offsets depicted in the top left and bottom left of Fig. 2, respectively, at 308W in the Atlantic Ocean.(top right) Neutral density g n at the same longitude.(bottom left) Thermodynamic neutral density based on the reference pressure depicted in the top-right panel of Fig. 2. (bottom right) McDougall and Jackett (2005b) materially conserved approximation to g n .

FIG. 2 .
FIG. 2. (top left) The latitude-depth dependent reference pressure seen by GPPD depicted in the top panel of Fig. 1. (top right) The reference pressure associated with Lorenz reference state underlying thermodynamic neutral density depicted in the bottom left of Fig. 1. (bottom left) Density offsets entering the construction of GPPD.(bottom right) Scatterplot of the GPPD reference pressure against GPPD density offset, showing the linear dependence of one on the other.

FIG. 3 .
FIG.3.The piecewise polynomial function of z R entering the construction of s r1d (z r ) (see text for details; thin solid line).Location of the maximum of the joint pdf between g n 2 s[S, u, p r (z r )] 1 ajz r j 1 b and jz R j (blue and orange dots).

FIG. 5
FIG. 5. (top left) Materially conserved part of g n obtained by bin-averaging g n in (u, S) space.(top right) The quasi-material Lorenz neutral density bin averaged in the same way as g n .(bottom left) Difference between (top left) and (top right).(bottom right) As in bottom left, but for g a .

FIG. 6 .
FIG. 6. Probability distribution functions for the sine of the angle between =g T and the neutral vector compared with that of various other density variables: (a) s 0 , (b) s 2 , (c) g a , and (d) g n .