Estimating Lorenz’s Reference State in an Ocean with a Nonlinear Equation of State for Seawater

ThestudyofthemechanicalenergybudgetoftheoceansusingtheLorenzavailablepotentialenergy(APE)theoryisbasedonknowledgeof theadiabaticallyrearrangedLorenzreferencestateofminimumpotentialenergy.Thecompressibleandnonlinearcharacteroftheequationofstateforseawaterhasbeenthoughttocausethereferencestatetobeilldeﬁned,castingdoubtontheusefulnessofAPEtheoryforin-vestigating ocean energetics under realistic conditions. Using a method based on the volume frequency distribution of parcels as a function of temperature and salinity in the context of the seawater Boussinesq approximation, which is illustrated using climatological data, the authors show that compressibility effects are in fact minor. The reference state can be regarded as a well-deﬁned one-dimensional function of depth, whichformsasurfaceintemperature,salinity,anddensityspacebetweenthesurfaceandthebottomoftheocean.Foraverysmallproportion ofwater masses, this surfacecan bemultivalued and water parcels can have upto twostatically stable levels in the referencedensity proﬁle, of which the shallowest is energetically more accessible. Classifying parcels from the surface to the bottom gives a different reference density proﬁle than classifying in the opposite direction. However, this difference is negligible. This study shows that the reference state obtained by standard sorting methods is equivalent to, though computationally more expensive than, the volume frequency distribution approach. The approachthatispresentedcanbeappliedsystematicallyandinacomputationallyefﬁcientmannertoinvestigatetheAPEbudgetoftheocean circulation using models or climatological data.


Introduction
The global kinetic energy (KE) budget plays a key role in ocean energetics, for it is often the natural starting point for discussing how the ocean circulation is Denotes Open Access content. forced and dissipated. In its standard form, the kinetic energy budget reveals that kinetic energy is primarily controlled by 1) the power input due to the wind forcing and tidal forcing, 2) viscous dissipation, and 3) the net conversion between potential energy (PE) and kinetic energy (Gregory and Tailleux 2011). Of these three terms, only the first one is well understood, and quantifiable from observations, as discussed in Roquet et al. (2011) and Zhai et al. (2012), among others. Although local turbulent viscous dissipation rates are routinely observed as part of campaigns to measure turbulent mixing in the oceans, these are so variable spatially and temporally that it is difficult to infer what the volume-integrated viscous dissipation is. Oort et al. (1994) originally assumed that in the ocean there is a net conversion of potential energy into kinetic energy. This was then disputed by Toggweiler and Samuels (1998) and in subsequent modeling studies using ocean-only models (Gnanadesikan et al. 2005) and coupled climate models (Gregory and Tailleux 2011), in which net conversion from kinetic energy to potential energy has been found, dominated by wind power input through Ekman pumping in the Antarctic circumpolar region.
The standard kinetic energy budget, however, does not directly reveal the role and importance of surface buoyancy fluxes and interior mixing processes of heat and salt in forcing and dissipating the ocean circulation. These processes control, at least in part, the sign and magnitude of the net conversion between potential and kinetic energy in the oceans. Linking the mechanical energy (gravitational potential plus kinetic energy) budget to surface buoyancy fluxes and interior mixing processes has been a controversial topic and has been overlooked in recent reviews on ocean energetics by Wunsch and Ferrari (2004), Kuhlbrodt et al. (2007), and Ferrari and Wunsch (2009). In these reviews, the mechanism by which surface buoyancy forcing by heat and freshwater fluxes provides mechanical energy to the ocean circulation has been related to the mechanical work of expansion and contraction produced by a heat engine, which in classical thermodynamics can be done by combining energy conservation with the entropy budget. For a stratified fluid, however, the available potential energy (APE) budget introduced by Lorenz (1955) is better suited for representing the mechanical power input by buoyancy forces that results from surface heat and freshwater fluxes. Also, available potential energy decomposes the net conversion between potential and kinetic energy as the difference between a production term by surface buoyancy fluxes minus a dissipation term by interior mixing (Hughes et al. 2009;Tailleux 2010Tailleux , 2012. Quantifying available potential energy in the framework of Lorenz (1955) relies on the definition of a globally static and stably stratified reference state obtained from the actual state by means of an adiabatic rearrangement of the fluid parcels conserving mass and salt. Defining the reference state for a fluid with a linear equation of state (EOS) in a simple domain with no sills or enclosed basins is made straightforward by, for instance, sorting water parcels according to density and filling the ocean volume level by level with the sorted density field (e.g., Winters et al. 1995). This sorting formulation has been used recently to investigate the mechanical energy budget of ocean circulation in a number of idealized situations (Hughes et al. 2009;Saenz et al. 2012;Hogg et al. 2013;Dijkstra et al. 2014).
Generalizing the concept of a reference state for the ocean under realistic conditions has remained problematic because of the presence of topographic sills and because of the nonlinearities of the equation of state for seawater. Topographic barriers block heavy waters from flowing between ocean basins, and it has been unclear how to represent this effect on the APE budget. Stewart et al. (2014) show that energy fluxes in the APE budget are largely insensitive to how the effects of topographic barriers are represented. On the other hand, the compressibility and the nonlinear dependence that density has on pressure leads to difficulties when sorting water parcels according to density, as the density needs to be recalculated at every level that is being filled. This makes calculating the reference state computationally expensive and causes the position of water parcels in the reference state to depend on the number of levels used and on whether the ocean volume is being sorted by increasing or decreasing density (e.g., Ilicak et al. 2012). Ad hoc sorting methods have nevertheless been devised to overcome these difficulties and have been used to investigate, for example, available potential energy in the global oceans (Huang 2005) and mixing in fluids with a nonlinear EOS (Ilicak et al. 2012;Petersen et al. 2015;Butler et al. 2013), but they provide no physical insight into the effects of compressibility and nonlinearities on the resulting reference density profile and remain computationally expensive. The lack of a physically tractable and well understood method to construct the Lorenz reference state has led to skepticism about the applicability of the APE framework to the ocean with a nonlinear equation of state. However, the physical basis for such skepticism has been lacking.
In this paper, we investigate the reference state for the ocean with a nonlinear equation of state for seawater by generalizing the approach proposed by Tseng and Ferziger (2001), based on the volume frequency distribution of water masses in temperaturesalinity space. We demonstrate that, for all practical purposes, the oceanic reference state can be regarded as well defined. In doing so, we provide a framework by which the reference state can be obtained and characterized systematically in a physically tractable and computationally efficient manner. The work by Stewart et al. (2014) and the work we present in this paper together lay the groundwork for using available potential energy to quantify the mechanical energy budget in ocean circulation under realistic conditions.
In section 2, we summarize basic concepts of the local formulation of APE theory for a Boussinesq ocean that will be used in this paper. Then, in section 3, we present a new approach based on the volume frequency distribution of water masses to link the reference state to the thermophysical properties of the ocean. In section 4, we discuss the reference position of the ocean parcels obtained using our approach, and in section 5, we discuss the links of our approach with sorting-based methods. Last, in section 6, we briefly discuss the implications for ocean APE estimates, and we summarize and discuss the results in section 7.

A review of available potential energy for a Boussinesq fluid
To set the context and to introduce concepts that will be used in this paper, in this section we present a brief summary of the general theoretical framework in Tailleux (2013) and present it in a form that is applicable to the ocean.

a. Model assumptions and governing equations
We investigate the Lorenz reference state for a realistic ocean in the context of the Boussinesq approximation. The governing equations are where u 5 (u, y) is the horizontal velocity field; w is the vertical velocity; D/Dt 5 ›/›t 1 u›/›x 1 y›/›y 1 w›/›z; $ h 5 (›/›x, ›/›y); p is the pressure; F u represents momentum forcing; r is density; r 0 is a reference Boussinesq density; g is the acceleration of gravity; f is the Coriolis parameter; Q is the potential temperature or Conservative Temperature (as in McDougall 2003); S is salinity; and _ Q and _ S represent sinks and sources for Q and S, respectively. The equation of state r 5 r[S, Q, p 0 (z)] is assumed to be formulated in terms of the Boussinesq pressure p 0 (z) 5 2r 0 gz to ensure energy consistency (e.g., Young 2010;Tailleux 2012), which also requires the surface to be situated at z 5 0.
b. Available potential energy: Work by buoyancy forces By construction, the reference state in APE theory is a notional state of rest that can, in principle, be reached from the actual state by means of a volume-conserving, adiabatic, and isohaline rearrangement. The reference density r r (z, t) is thus a one-dimensional function of depth and is also a function of time in the sense that the ocean can have different reference states at different times. It is in hydrostatic equilibrium with the reference pressure p r , that is, ›p r /›z 5 2r r g. As discussed in Tailleux (2013), it is natural in APE theory to define the buoyancy of a fluid parcel relative to the reference density r r , that is, The APE for a fluid parcel E a can then be defined as the work against buoyancy forces required to lift a parcel from its equilibrium position in the reference state z r to its actual position z. We refer to E a as APE density, after Winters and Barkan (2013). For a fully compressible fluid (Andrews 1981) and for a Boussinesq fluid with a linear equation of state (Holliday and Mcintyre 1981), the APE density can be written as By construction, the reference position z r of a fluid parcel in the reference state corresponds to the depth at which its density is equal to the reference density: which is referred to as the level of neutral buoyancy (LNB) equation. Equation (8) defines z r 5 z r (S, Q, t) as a function of materially conserved variables and time.
For small-amplitude displacements z 5 z 2 z r , E a reduces to the familiar formula where N 2 r is the squared buoyancy frequency: in which c 2 s 5 (›r/›p) 21 is the squared speed of sound. Note that while r r (z, t) depends only on depth for any given time t, N r (x, t) 5 N r [S(x, t), Q(x, t), t] depends on position and time. The quantity N 2 r will be useful later in this paper to select a root when the LNB equation has more than one solution.

c. Interpretation of APE density and its relation to global APE
Integration by parts of the first term in the right-hand side of Eq. (7) allows the change in potential energy following a parcel between its actual and reference position to be written as ð z z r gr[S, Q, p 0 (z 0 )] dz 0 5 rgz 2 r r gz r 1 where gpe 5 rgz/r 0 and ie are specific gravitational potential and internal energies in the Boussinesq approximation, respectively. In a fully compressible fluid, the change in specific internal energy associated with a change in pressure is given by die 5 2pdn adiab 5 pdp/(r 2 c 2 s ), where n 5 1/r is specific volume. For a Boussinesq fluid, this change is approximated as die bouss ' (2r 0 gz)(2r 0 gdz)/(r 2 0 c 2 s ) 5 g 2 zdz/c 2 s .
Thus, the term involving c 2 s in Eq. (7) is interpreted as the change in specific internal energy Die between the parcel's actual position and its position in the reference state. As a result, and since the reference state is in hydrostatic equilibrium, the APE density can be written as r 0 E a 5 r 0 (Dgpe 1 Die) 1 p r (z, t) 2 p r (z r , t) , (12) that is, as the sum of changes in gravitational potential and internal energy, as well as pressure work, as originally established by Andrews (1981) for a fully compressible fluid. The volume integral of p r (z, t) 2 p r (z r , t) vanishes for reference states that are rearranged in an adiabatic and isohaline manner from the actual state (see, e.g., Tailleux 2013). As a result, the volume-integrated APE density becomes ð that is, the difference between the potential energy of the actual and reference states as defined by Lorenz (1955). This establishes the equivalence between the local approach presented in this section and the global approach used by Lorenz (1955). Note that, from the above result, total APE can also be written as the sum of its gravitational potential and internal energy components, exactly as for a fully compressible fluid. The internal energy contribution arises because of the pressure dependence of density, which explains why it is absent in idealized studies using a linear equation of state (e.g., Winters et al. 1995;Hughes et al. 2009). Further, the Boussinesq approximation for seawater admits a fully consistent set of thermodynamic functions that can be obtained from the exact Gibbs function describing all thermodynamic properties (for details see Tailleux 2012).

d. Local description of the ocean energy cycle
An advantage of the APE density in Eq. (7) over global APE is that it gives a local description of KE/PE conversions. The local evolution equation for the kinetic energy is obtained by multiplying the horizontal momentum equations by u (after subtracting p r from p), which leads to where E k 5 u 2 /2 is the kinetic energy density of the circulation (which neglects the vertical velocity component in the hydrostatic approximation). A local evolution equation for the APE density is obtained by taking the time derivative of Eq. (7) premultiplied by r 0 . After using the definition of the LNB and some manipulation (for details see Tailleux 2013), we obtain (16) In the above equation, G a represents the local generation of APE, given by where G Q and G S are thermal and haline efficiency factors defined by where a and b are the thermal expansion and haline contraction coefficients defined relative to the (Q, S, p) variables, and a and b are their vertically averaged values over the interval [z r , z]. Alternatively, using u Á $ h (p 2 p r ) 5 $ Á [(p 2 p r )v] 1 w(r 2 r r )g, the local energy cycle can also be written as Although the two formulations are mathematically equivalent, they provide alternative expressions for the local conversion between APE and KE, C(A, K). In the first formulation, this conversion is given by while in the second formulation it is expressed in terms of the buoyancy flux: So far, most studies of ocean energetics have used Eq. (23) as the definition of C(A, K), the exception being Gregory and Tailleux (2011), who used Eq. (22) to investigate the spatial and temporal variations of the vertically integrated C(A, K) in global warming experiments.

e. APE production and dissipation
A common approach in ocean modeling is to express the source and sink terms for temperature and salinity as the divergence of heat and salt fluxes, namely, _ Q 5 $ Á F Q and _ S 5 $ Á F S , associated with the surface boundary conditions F Q Á n 5 Q/(r 0 c p ) and F S Á n 5 S 0 (E 2 P), where Q is the surface heat flux into the ocean (with Q . 0 denoting heat going into the ocean), c p is the heat capacity, S 0 is a reference salinity value, and E 2 P is evaporation minus precipitation. Assuming insulating boundaries everywhere else, the volume integral of G a in Eq. (17) can be decomposed as the difference between a production term by surface buoyancy fluxes G(A) minus a term representing dissipation by diffusive fluxes of heat and salt D( The APE production term is given by whereã andb are the mean values of a and b applicable to each parcel at the surface. Since z r , 0 by definition, it follows that APE production is associated with local cooling (Q , 0) and net evaporation (E 2 P , 0). For a linear equation of state and in the absence of salinity variations, we have G(A) 5 2 Ð S gaz r Q/c p dS, which reduces to the term denoted by F b2 in Hughes et al. (2009). Similarly, the dissipation term is 3. Linking the reference state r r (z) to the stratification and thermophysical properties of the ocean In this section, we outline a new approach to constructing the Lorenz reference state r r (z). This approach can be regarded as an extension to a binary fluid of the approach previously developed by Tseng and Ferziger (2001). As a result, it is expected to be computationally more efficient than the commonly used sorting methods (e.g., Winters et al. 1995;Huang 2005). We begin by making the important point that in practice determining r r (z) only requires the knowledge of r r (z k ) at a relatively small number N of target depths z k , k 5 1 . . . N. Typically 30 # N # 50 is sufficient to accurately describe the ocean stratification. Assuming that we can devise a computationally inexpensive way to estimate r r (z k ) at the preselected N target depths z k (to be discussed below), simple interpolation can then be used to estimate r r (z) at any arbitrary depth z if needed. Once r r (z) is known, the reference position z r (x) for each ocean parcel can be obtained by solving the LNB Eq. (8), for instance using a standard numerical root-finding scheme. In fact, the LNB Eq. (8) defines z r 5 z r (Q, S) as a function of the materially conserved variables Q and S, so that we can From a fundamental viewpoint, the above method divorces the task of constructing the reference density profile r r (z) from the task of computing the reference position z r (x) for each parcel of water. This is an important departure from current sorting-based approaches, which accomplish the two tasks simultaneously. As we justify below, such a divorce is not only feasible but also essential for elucidating the nature of the difficulties associated with the nonlinearity of the equation of state and its binary character through the dependence on temperature and salinity. From a computational viewpoint, such a divorce also proves essential for understanding how to design algorithms that are computationally more efficient and parallelizable than sorting-based approaches.
In this section, we describe the construction of the reference density profile. We start in section 3a by deriving a general equation of the form f [r r (z)] 5 0, which can be easily solved numerically for the reference density profile r r (z). In section 3b, we discuss how compressibility may lead to a special situation that may arise when solving f [r r (z)] 5 0 for r r (z). We then describe the two choices that can be made for arriving at a solution of this equation in section 3c.
a. An implicit equation for r r (z) based on linking volume in physical space and thermohaline space We now set out to define r r (z) by linking the volumetric properties of the ocean stratification to the thermophysical properties of seawater in temperaturesalinity (thermohaline) space. This is done for a Boussinesq ocean at some instant in time t, with potential or Conservative Temperature Q(x, y, z) and practical or Absolute Salinity S(x, y, z) assumed to be given as functions of longitude x, latitude y, and depth z increasing upward, with the ocean surface being located at z 5 0. For brevity, we omit the time dependence in the following. In this paper, the density of a fluid parcel is estimated using the nonlinear equation of state of Jackett and McDougall (1995), which is defined in terms of potential temperature, practical salinity, and pressure.
As we illustrate below, the idea underlying our approach is that the volume of seawater V(z l , z u ) occupied between two arbitrary levels z l # z u can be written equivalently either in physical or thermohaline space. In physical space, V(z l , z u ) can be written as where A(z) is the area of the ocean at a given depth z, which can be obtained from the knowledge of the ocean bathymetry by defining depth H(x, y) as a function of horizontal position (x, y). To understand how to obtain an expression for V(z l , z u ) in thermohaline space, it is necessary to construct the normalized volume frequency distribution function p(Q, S), such that p(Q 0 , S 0 )dQdS represents the volume fraction of the water masses in the region of the temperature-salinity domain [S 0 , S 0 1 dS] 3 [Q 0 , Q 0 1 dQ] and satisfies the normalization condition where the thermohaline domain [S min , S max ] 3 [Q min , Q max ] is assumed to contain all possible Q, S values encountered in the ocean. In other words, the normalized volume frequency distribution function p(Q, S) is constructed by adding up the volume of all grid cell boxes of a discrete global ocean grid with values of Q, S within the range [S 0 , S 0 1 dS] 3 [Q 0 , Q 0 1 dQ] and then dividing by the total ocean volume. Assuming that r r (z u ) is known, then V(z l , z u ) must by construction occupy the volume between the two isopycnals of equations r[S, Q, p 0 (z l )] 5 r r (z l ) and r[S, Q, p 0 (z u )] 5 r r (z u ) in (Q, S) space. In thermohaline space, an isopycnal r at constant pressure p is best viewed as the locus in (Q, S) space satisfying the equation S 5 S(Q, r, p), where the function S(Q, r, p) can be obtained by inverting the equation of state r 5 r(S, Q, p), as illustrated in Fig. 1. Note that because of the nonlinearity in the function r 5 r(S, Q, p), a well-defined, one-to-one mapping function Q 5 Q(S, r 0 , p) cannot be obtained. In this paper, the function of state S 5 S(Q, r, p) is referred to as a ''salinity curve'' and is obtained by inverting the equation of state by Jackett and McDougall (1995). This makes it possible to describe the isopycnals r 5 r r (z l ) and r 5 r r (z u ) as two salinity curves with equations S 5 S[Q, r r (z l ), p 0 (z l )] 5S(Q, z l ) and S 5 S[Q, r r (z u ), p 0 (z u )] 5S(Q, z u ), respectively. Note that because both r r (z) and p 0 (z) depend on z only, the salinity curve of equation S 5 S[Q, r r (z), p 0 (z)] 5S(Q, z) can be viewed as a function of Q and z alone, which is indicated by the tilde over S. As a result, a mathematical expression for the volume V(z l , z u ) in thermohaline space can be written as where V T is the total ocean volume, defined by while the function s(Q; z l , z u ) represents the volume fraction contained between z u and z l at fixed Q and is given by The justification for Eq. (31) arises from the fact that while the two salinity curvesS(Q, z l ) andS(Q, z u ) are in general well separated, and such thatS(Q, z l ) . S(Q, z u ) whenever z l , z u , these two curves may nevertheless occasionally intersect in regions where the reference level z r (Q, S) is multivalued and associated with multiple LNBs, as we discuss later. When the two salinity curves intersect, their intersection point (Q c , S c ) is singular and can be shown to be associated with a vanishing squared buoyancy frequency N 2 r (Q c , S c ) 5 0. Because the integral Eq. (29) becomes ambiguous in regions where salinity curves intersect (discussed in section 3b), defining s as per Eq. (31) is one way to avoid the ambiguity and the double counting of water masses and allows the bottom-up and top-down approaches discussed below to be connected with the bottom-up and top-down sorting strategies discussed later in the text.
Equation (29) is a key result of this paper, for it shows that if r r (z u ) is known, then Eqs. (29) and (27) together provide an implicit equation of the form which can be easily solved for r r (z l ) using a numerical root searching method. The overall approach is feasible because the reference density r r (z) is known a priori both at the ocean surface (z 5 0), where it is equal to the lowest possible density r min , as well as at the maximum ocean depth z 5 2H max , where it must be equal to the maximum possible density value r max . Therefore, the corresponding bounding salinity curves S 5 S[Q, r min , p 0 (0)] and S 5 S[Q, r max , p 0 (2H max )] in (Q, S) space are also known, and one may start at the surface by setting r r (z u ) 5 r min and integrate toward the bottom of the ocean. We construct a normalized volume frequency distribution function p(Q, S) using the annually averaged temperature and salinity fields from the 2009 version of the Levitus World Ocean Atlas (WOA09; Locarnini et al. 2010;Antonov et al. 2010). Data from the WOA09 were classified on a discretized Q, S plane with bin sizes dQ 5 0.0058C and dS 5 0.002 practical salinity units (psu). The WOA09 dataset is defined on a vertical grid consisting of 33 discrete levels separated by dz varying from 10 m at the surface to 500 m at the bottom of the ocean. We use the WOA09 vertical grid but with higher resolution in the first 10 m, where we defined levels every 1 m, resulting in a vertical grid with 42 levels. Values for thermodynamic constants were obtained from McDougall and Barker (2011), and we used gravity g 5 9.81 m s 22 and reference density r 0 5 1027.0 kg m 23 . The left-hand side of Eq. (29) was calculated using topography from the gridded global relief model of Earth's surface in Amante and Eakins (2009), assuming a sphere radius of R 5 6.4 3 10 6 m. We solve Eq. (29) at each of the 42 levels of our vertical grid by integrating from the surface toward the bottom of the ocean. The reference density profile obtained this way is shown in Fig. 2. In the remainder of this paper, we will discuss the physical properties of this profile.

b. Effects of compressibility: Intersection of salinity curves
For a given target depth z 5 z u , let us consider the isopycnal r 5 r r (z). As mentioned above, the latter is most conveniently described in terms of the salinity curve S 5 S[Q, r r (z), p 0 (z)]. Let us now move toward the next target depth obtained as described above, z 1 dz 5 z l , associated with the isopycnal r 5 r r (z 1 dz) and salinity curve S 5 S[Q, r r (z 1 dz), p 0 (z 1 dz)]. We are now interested in what controls the separation between these two salinity curves and, in particular, whether these two curves can intersect and what the conditions are for them to intersect.
Using the approximations r r (z 1 dz) ' r r (z) 1 [›r r (z)/›z]dz and p 0 (z 1 dz) ' p 0 (z) 2 r 0 gdz, it follows that which implies that, to leading order, the distance dS 5S(Q, z 1 dz) 2S(Q, z) is given by where c 2 s is the squared speed of sound. As a result, it follows that where N 2 r , the squared buoyancy frequency, N 2 r 5 2(g/r 0 )[(›r r /›z) 1 r 0 g/c 2 s ] [Eq. (10)]. Equation (36) indicates that the distance separating two isopycnals and the tilt between them are controlled by the stability of the reference state. Since c 2 s and N 2 r are functions of Q, S and p 0 (z), they can vary along a given isopycnal, as can dS, causing the isopycnal r r (z 1 dz) to tilt with respect to isopycnal r r (z). For a given tilt, if the separation between the two isopycnals is large enough, the isopycnals r r (z) and r r (z 1 dz) are well separated in the sense that they do not intersect, as shown in Fig. 1. However, if the   FIG. 2. (top) Reference density profile r r obtained with the top-to-bottom sorting scheme using the highest resolution possible. Also shown are the reference density profiles obtained with the volume frequency distribution method using the vertical resolution of the WOA09 data (with resolution varying between Dz 5 1 m at the surface and 500 m at depth), 10 and 0.1 m as indicated. (bottom) Relative error with respect to the top-to-bottom sorting scheme using the highest resolution possible: 100[(r r 2 r r,sorted )/r r,sorted ]. The relative error is 0 at the surface and at the bottom (by construction). Note that the limits of the horizontal axis are set to show the range of values of relative errors that occur below z 5 220 m. Within 20 m of the surface, the relative error peaks to about 0.45% for 42 levels upward and downward, 6.5 3 10 23 % for dz 5 10 m, and 0.75% for dz 5 0.1 m.

MAY 2015 S A E N Z E T A L .
tilt is large enough, the isopycnals could intersect as illustrated in the case shown in Fig. 3. The separation between isopycnals is dictated by Eq. (29) and by the relative volume of the ocean contained between r r (z) and r r (z 1 dz) as illustrated in Fig. 4. For the WOA09 annual-mean climatology, Fig. 4 shows that a large volume of the ocean is concentrated in a very small region in Q, S space. When this large volume of the ocean is somewhere between isopycnals r r (z) and r r (z 1 dz) in Fig. 3, the separation between these isopycnals, specifically dz, does not have to be very large for Eq. (36) to be satisfied. Mathematically, N 2 r (Q, S) 5 0 defines a caustic line in thermohaline space that separates the region N 2 r . 0, where the curves do not intersect, from a region where z r (Q, S) can have up to three possible values (this will be discussed in more detail in section 4).

c. Top-down versus bottom-up approaches
There are two ways to construct r r (z) based on Eqs.
(27) and (29). Since r r (z) is known at the surface, one may use a top-down approach by starting at the ocean surface and integrating downward. Alternatively, and since r r (z) is also known at the bottom of the ocean, a bottom-up approach may be implemented in which integration is carried out starting at the bottom and integrating upward. We will see later in section 5 that differences between these two approaches are small.

d. Topographic barriers
Topographic enclosures may trap dense fluids and isolate them from flowing into other ocean basins. The effects that topographic barriers may have on APE or the reference state are ignored here in the sense that all parcels are treated as if they were all in the same ocean basin. The reader is referred to Stewart et al. (2014) for a more detailed discussion of the effects of topographic barriers.

Position of a water parcel in the reference state
The reference position occupied by a fluid parcel is the depth z r at which the in situ density of the parcel brought adiabatically (and conserving salinity) to that level is equal to the reference density or r[Q(x, y, z), S(x, y, z), p 0 (z r )] 5 r r (z r ) .
From this relation it is clear that the reference position of the fluid parcel is independent of the path taken to reach that position and hence that the buoyancy forces are conservative. The depth level z r is referred to as the LNB by Tailleux (2013). Equation (37) defines z r (S, Q) as a function of S and Q and is illustrated graphically in Fig. 5. Two distinct regions can be identified in the reference state profile r r (z) obtained with the WOA09 annualmean climatology (Fig. 5, lower panel). In the upper 1 km of the reference state profile, the adiabatic change in in situ density with depth of a water parcel occupying a given level is significantly smaller than the rate at which the reference density increases with depth: j›r/›zj 5 jr 0 g/c 2 s j ( jdr r /dzj, so that N 2 r ' 2(g/r 0 )r 0 r (z) is dominated by the vertical derivative of r r (z). According to Eq. (36), this corresponds to a regime where the salinity curves and isopycnals are well separated. In the region where 24 , z , 21 km on the other hand, j›r/›zj Q,S j 5 r 0 g/c 2 s becomes of the same order as jr 0 r (z)j. As a result, N 2 r becomes relatively small and can even vanish, resulting in the salinity curves and isopycnals approaching each other and even intersecting, as illustrated in Figs. 3 and 4.
The condition N 2 r 5 0 defines a caustic line in (Q, S), represented by the black line in Fig. 6, which separates the thermohaline space into two regions, the first where the LNB Eq. (37) has only one possible stable solution with N 2 r . 0 everywhere and the second where multiple solutions to the LNB equation exist. In the latter region (which we will refer to as the overlap region, seen in the upper panel of Fig. 5 and detailed in Fig. 6), we find that in general the LNB equation has up to three possible solutions, two being stable and associated with N 2 r . 0, and one unstable and associated with N 2 r , 0. More specifically, we find a stable and shallow solution z r1 , an unstable intermediate solution z r2 , and a stable deep solution z r3 . On the edges of the overlap region, Eq. (37) has two roots z r1 and z r2 , only one of which is statically stable, namely, the deepest one z r2 .
The existence of the overlap region and two stable solutions explains why, in general, the bottom-up and top-down algorithms discussed previously yield two distinct reference density profiles. This is because the top-down approach naturally selects the stable solution with the largest z r (which tends to have larger N 2 r ). As a result, the top-down approach can be regarded as yielding a more stable reference state than the bottomup approach.
Based on our analysis of the WOA09 annual-mean climatology, we find that the volume of water masses found in the overlap region represents a negligible fraction of the overall ocean volume. It follows that the reference position for a large majority of the fluid parcels is unique and stable. The volume of the water parcels inside the overlap region is 0.0023% of the global ocean volume in the WOA09 dataset. The vast majority of these water parcels are Antarctic surface waters, in the Southern Ocean south of 608S, with the remaining volume located in the Arctic Ocean (Figs. 7, 8). All parcels with three possible reference levels are within the upper 600 m of the surface of the ocean, most of them above 200 m (Fig. 7). In other words, their actual location is at much shallower depths than their shallow reference level z r1 between the depths of 1.5 and 2.5 km (Figs. 5, 6). For these water parcels to reach the deep reference level z r3 , they would need to use extra energy to go below z r2 . Therefore, the shallow level is more energetically accessible, and we choose z r 5 z r1 as the reference level for the parcels in the overlap region.
Cabbeling, the densification that occurs when water masses with different Q, S mix, results from the fact that the isodensity curves in Q, S space are curved. Mixing of any pair of water parcels that lie on the same isopycnal can lead to some level of cabbeling, regardless of where these parcels lie on the reference state curve. Our analysis gives an indication of the conditions under which densification from cabbeling can have the largest effects on local and global APE and on the dynamics of the circulation in the ocean. Cabbeling from mixing between water masses that have a position on the reference state curve inside the overlap region and water masses of slightly higher salinity, with a position on the reference state curve outside of this region (Fig. 6), can lead to drastic changes in z r . Large changes in z r resulting from the mixing of these parcels would lead to large changes in buoyancy [Eq. (6)] and APE density [Eq. (7)] in the regions of the ocean shown in Figs. 7 and 8.
By systematically implementing the procedure to obtain the reference level described so far, we are able to calculate z r (x, y, z) for the WOA09 annual-mean climatology. The results of these calculations have a remarkable resemblance to the calculations by Saenz et al. (2012) for an idealized ocean basin with realistic forcing. Water parcels with the deepest reference level (lowest z r ) and the largest available potential energy density E a defined in Eq. (7) are at the surface of the Southern Ocean in the Weddell Sea close to the Antarctic Peninsula (left and right upper panels of Figs. 9 and 10, respectively) and in the western part of the Ross Sea (left and right upper panels of Fig. 9 and lower panels of Fig. 10, respectively), where z r , 24 km. In the rest of the Southern Ocean, the reference levels vary between 21 and 24 km, with most of the surface waters having z r toward the shallower end of this range of depths. This is consistent with the regions where Antarctic Bottom Water (AABW) is generated and with the depths at which this water mass is found in the ocean. In the Arctic Ocean, the regions with the deepest reference level and highest available potential energy density are in the Greenland Sea and the Barents Sea, where z r can be deeper than about 22.5 km, consistent with the regions where North Atlantic Deep Water (NADW) originates and with the depth at which this water mass is found in the ocean. The region of confluence between the Arctic Ocean and the North Atlantic Ocean and the Labrador Sea (Figs. 9, 10) has intermediate values of available potential energy density and reference levels of up to about 21 km. In the Mediterranean Sea, off the coast of France and in the Adriatic Sea, the reference level can reach depths of 21 and 21.5 km, respectively. In the rest of the surface of the ocean, the reference level is within 500 m of the surface of the ocean (Fig. 10). The spatial distribution of z r is broadly consistent with the spatial distribution of available potential energy density E a .
The rate of APE production by surface buoyancy fluxes is given by Eq. (25): which reduces to the term F b2 in Hughes et al. (2009) for a linear equation of state. Our results suggest that in the regions shown in Fig. 9, where water parcels have a large z 2 z r and where surface buoyancy fluxes are large (i.e., in the Weddell Sea close to the Antarctic Peninsula, in the western part of the Ross Sea, and in the Greenland Sea and Barents Sea), there can be large sources of available potential energy, and thus of mechanical energy, to the global overturning circulation (Saenz et al. 2012).

The reference configuration obtained by sorting water parcels
The procedures to calculate r r (z) and z r discussed so far can also be understood by interpreting them from the perspective of adiabatic and isohaline sorting of water parcels. Consider a scheme that sorts water parcels, where each water parcel corresponds to a grid cell volume from a discretized ocean volume, and assume that sorting is done from the top of the ocean to the bottom, assigning the lightest available water parcel to each level. Consider the situation in which the water column has been sorted from the surface down to cell k at level z 5 z k , where it has density r k . The pressure on the upper face of level k 1 1 is p 5 2gr 0 (z k 1 0:5Dz k ), where Dz k is the thickness of level k after being ''flattened'' out to occupy the volume of the ocean at that level. This pressure is used to calculate the density that the unsorted parcels would have if placed at level k 1 1. The water parcels that remain to be sorted are those for which S # S[Q, r r (z k ), z k ] in Fig. 1. The parcel with the smallest density, parcel m, is assigned to level k 1 1. The temperature and salinity of this parcel Q m , S m correspond to a class with volume fraction p m 5 p(Q m , S m ) 5 V m /V T , assuming that each water parcel has a unique p(Q, S) class [this is true for the WOA09 annual-mean climatology we analyzed; however, parcels with the same p(Q, S) class can be lumped together into a single parcel]. Assigning the parcel with the smallest density to level k 1 1 is equivalent to satisfying the discrete version of Eq. (29) with the most accurate vertical step size Dz possible, dictated by the volume of grid cell m divided by the area of the ocean at level k, A(z k ).
The special case (extremely unlikely, however, as we have found with the WOA09 data) in which more than one water parcel has the same density at level k 1 1 is dealt with as follows. Assume that there are several water parcels with different Q, S values that have the same density at level k 1 1. On a given isopycnal, compressibility increases monotonically with decreasing Conservative Temperature [as illustrated in Fig. 11, calculated using subroutines in McDougall and Barker (2011)]. Sorting these water parcels by compressibility and placing the least compressible one at the shallowest level and the most compressible one at the deepest level results in the lowest potential energy state possible.
By assigning a parcel to a given level k, an adiabatic sorting scheme is equivalent to solving Eq. (37) for r r (z r ) and obtaining z r simultaneously. When a water parcel in the overlap region is chosen to occupy a level, it cannot be included in any of the remaining sorting iterations. In terms of the continuous integration of Eq. (29), this means that when solving for r r (z l ), regions with dS , 0 inside the overlap region (Fig. 6) should not be included in the integration. Thus, when sorting from the surface toward the bottom of the ocean, water parcels in the overlap region are assigned to their shallowest reference level. Similarly, in solving Eq. (29) by integrating from the bottom of the ocean toward the surface, steps with dS . 0 should not be taken. Sorting in this direction assigns water parcels with Q, S in the overlap region to their deepest reference level. Recall that the shallowest stable reference level is the most energetically accessible. Therefore, sorting from top to bottom is preferred. As a result, when the above restrictions are enforced, the reference density curve obtained by integrating Eq. (29) with dz . 0 results in a profile different to when the integration is performed with dz , 0 (Fig. 2). However, since the volume of the ocean within the overlap region is so small, this difference is negligible, and the profiles are indistinguishable. The potential energy of the reference state PE r in Eq. (13), or background potential energy, obtained for the WOA09 annual-mean climatology by sorting from the surface toward the bottom is 28.394 82 3 10 25 J, which is 1.4 3 10 24 % higher than PE r obtained by sorting from the bottom toward the top. This difference is well within the uncertainty of the data of any available dataset, either from measurements or simulations.
If the top-to-bottom sorting scheme and the volume integration procedure solving Eq. (29) each use the highest possible vertical resolution, given by the volume of each individual water parcel grid cell divided by the area of the ocean at each level, the difference between the two methods should be accounted for by the differences induced by not including parcels in the overlap region more than once in the sorting scheme. This is confirmed in Fig. 2, which shows the reference density profile obtained using the sorting scheme at the highest resolution possible, along with the four other reference density profiles obtained as follows: with the volume frequency distribution method by integrating from the bottom toward the top onto 42 levels and with the volume frequency distribution method by integrating from the top toward the bottom onto 42 levels and onto levels spaced at 10 and 0.1 m. Since the overlap region obtained using the WOA09 data is so small, the difference between the reference density profiles obtained with high vertical resolution with each method is very small, as shown in the bottom panel of Fig. 2. The error induced by coarser-resolution calculations is also small.
Existing sorting schemes are computationally very expensive. For example, the most accurate sortingbased approach one can conceive of uses as many target depths N p as there are parcels to sort, namely, the total number of grid cells in the discretized ocean. Thus, the algorithm involves N p 2 1 iterations at which the number of candidate points at the ith iteration is N p 1 1 2 i and involves O[N p (N p 1 1)/2] calls to the equation of state. Fewer target depths can be used to overcome this computational cost, at the expense of uncertainty and accuracy. This is to be contrasted with the two-step approach we present here. Using the volume integration procedure to obtain r r (z), solving Eq. (29) for a level k 1 1 requires the EOS to be evaluated only for the parcels within a narrow region in Q, S adjacent to the isopycnal curve at level k. Once r r (z) is known, only O(10) calls to the equation of state or less are required to solve the LNB equation [Eq. (8)] for each parcel, and MN p ( N p (N p 1 1)/2 calls to the equation of state are needed overall in the second approach. In addition, it is straightforward to implement the procedure of solving the LNB for each parcel in parallel within any existing ocean general circulation model. As a result, the volume integration procedure of solving Eq. (29) is a much faster computational procedure.

Consequence for estimates of the ocean APE
By considering a fully compressible ocean, including possible volume changes upon adiabatic rearrangement, Huang (2005) previously estimated the oceanic average APE to be 624.2 J m 23 . Since the total ocean volume is approximately V ocean ' 1:3 3 10 18 m 3 , this implies that the total APE is 8.11 3 10 20 J. In contrast, by computing the APE density and estimating its volume integral, we estimate that the total APE is 1.08 3 10 21 J, which yields an average APE value of 834 J m 23 (assuming the same value for V ocean ). The estimate by Huang (2005) is about 25% lower. The reason for this is unclear, however, FIG. 11. Conservative Temperature and Absolute Salinity diagram for in situ density (kg m 23 ; contours) and isentropic compressibility (Pa 21 ; color) at a pressure of 2000 db, from the International Thermodynamic Equation of Seawater-2010 (TEOS-10).

MAY 2015
S A E N Z E T A L .
because our estimate is also based on the Levitus annual-mean climatology (although for the WOA09 database). From a numerical viewpoint, the estimate for the total APE in Huang (2005) is based on computing the difference in potential energy between the actual state and the reference state, that is, between two extremely large numbers, which is potentially ill conditioned and prone to numerical errors [see also the relevant discussion by Winters and Barkan (2013)]. In contrast, our approach is based on integrating only relatively small positive numbers and is therefore expected to be better conditioned numerically. Whether this is sufficient to account for the differences between our estimates for APE and that of Huang (2005) is unclear and warrants further study. This is beyond the scope of this work, as the computational issues associated with deriving accurate estimates of APE require a full separate study of a more technical nature.

Summary and conclusions
We have discussed the effects that compressibility and nonlinearities of the equation of state of seawater have on the reference state and have illustrated these effects using annual climatological data for temperature and salinity in the ocean. Most of the water parcels in the ocean have a single, well-defined reference level in the reference state. Variations in compressibility with temperature and salinity cause a very small volume fraction of the ocean to have up to two stable positions in the reference state profile. We argue that because these volumes are located at high latitudes and shallow depths, one reference level, the shallowest, is energetically more accessible.
Our formulation, through Eqs. (29) and (32), allows us to strongly constrain (to a desired level of precision) the range of temperature and salinity properties that a water parcel at a given level must have in the reference state. Because we know either the overall minimum (surface) or maximum (bottom) density in the reference state ocean, we can proceed to efficiently construct the reference density profile by progressively filling the ocean volume adiabatically either downward from the surface or upward from the bottom.
We show that the adiabatic sorting schemes are equivalent to classifying water masses using the volume frequency distribution in temperature and salinity space, provided that the latter accounts for the water parcels that have been assigned a location in the reference profile. Because of the multiple equilibrium positions that some water parcels may have, different methods can yield different reference states. The differences between these profiles are so small that, within the uncertainty of available data from simulations and from climatologies, the reference state profile can be thought of as well defined for all practical purposes. Uniqueness can be enforced if one chooses to define the reference state as that obtained by filling the ocean volume downward from the surface, which we have argued is more energetically accessible than the one resulting from filling the ocean upward from the bottom of the ocean.
Our work here, along with the results in Stewart et al. (2014), provides the foundation for the investigation of the available potential energy budget of the ocean circulation under realistic conditions.