## Abstract

The vertically integrated potential energy of an incompressible stratified fluid formulated in density coordinates can be simply written as a weighted vertical sum of the squares of the vertical displacements of density surfaces, a general expression valid for arbitrary displacements. The sum of this form of potential energy and kinetic energy is then a conserved quantity for the multilayer shallow water model. The formulation in density coordinates is a natural one to find the Lorenz reference state of available potential energy (APE). We describe the method to compute the APE of an ocean state and provide two applications. The first is the classical double-gyre, wind-driven circulation simulated by a shallow water model at high resolution. We show that the eddy kinetic and eddy potential energies are localized in regions of large gradients of mean APE. These large gradients surround an APE minimum found between the two gyres. The second is the time-mean World Ocean Circulation reconstructed from hydrography (*World Ocean Atlas*) and reference velocities at 1000 db from the Argo float program to obtain an absolute circulation. The total available potential energy exceeds the total mean kinetic energy of the World Ocean by three orders of magnitude, pointing out the very small Burger number of the circulation. The Gulf Stream, the Kuroshio, the Agulhas retroflection, and the confluence regions are four examples that confirm the shallow water model results that large gradients of mean available potential energy can be used as predictors for the presence of high eddy kinetic energy (obtained here from satellite altimetry).

## 1. Introduction

Available potential energy (APE) has become a central element in oceanography because the presence of mesoscale eddies in the ocean is caused by baroclinic instability of the large-scale flow, releasing the APE associated with the large-scale circulation (Gill et al. 1974; Stammer 1997; Chelton et al. 2007; von Storch et al. 2012; among many other studies). The discussion of ocean energetics initiated by Munk and Wunsch (1998) and Wunsch and Ferrari (2004) has emphasized the kinetic energy budget, but Hughes et al. (2009) and Tailleux (2010, 2013) have stressed the importance of the available potential energy budget driven by the surface buoyancy fluxes to quantify the conversion between kinetic and potential energy. A quasigeostrophic formulation of potential energy density quadratic in the *small* vertical displacements was used by Oort et al. (1989) to estimate the World Ocean’s APE. However, the quasigeostrophic assumption of small vertical displacements is not correct since vertical excursions of the density surfaces lying at depths larger than 2000 m in the interior can outcrop near the Antarctic continent. Strong mesoscale eddies, such as Gulf Stream rings, are also associated with 800-m undulations of isotherms (Richardson et al. 1979). The small displacement approximation also appears invalid for nonlinear internal waves (Scotti et al. 2006). Following the importance given to ocean energetics, new methods to compute more accurately available potential energy in the ocean have emerged using pressure or height coordinates. What is required is the potential energy available to drive motions, the difference between the actual potential energy and a reference state chosen by Lorenz (1955), to be a state of rest with flat-density surfaces obtained by an adiabatic (conservative) redistribution of the actual state. Huang (2005) adapted the sorting algorithm introduced by Winters et al. (1995) to compute it. Saenz et al. (2015) improved the determination of the reference state when the nonlinear equation of state is considered. Another concept, the dynamical potential energy, defined as the horizontal anomaly of potential energy, was also introduced by Roquet (2013) for the general circulation and by Klymak and Moum (2003) and Moum et al. (2007) for internal waves.

At the difference of these studies, the present paper makes use of density coordinates to calculate APE. The use of a Lagrangian vertical coordinate is particularly attractive because fluid parcels can be followed adiabatically between the reference state of rest and the observed state. The top-to-bottom integrated APE reduces to a weighted sum of quadratic vertical displacements of density interfaces from their Lorenz reference positions, which is valid for *arbitrary* vertical displacements. This form of potential energy was introduced by Ripa (1991) in the context of the derivation of stability theorems for the multilayer shallow water model formulation, but the value of this form of potential energy has been largely overlooked. Indeed, Hughes et al. (2009), von Storch et al. (2012), Zemskova et al. (2015), and others discuss the oceanic energy cycle in *z* or pressure coordinates. Isopycnal models [such as the Miami Isopycnal Coordinate Ocean Model (MICOM)], introduced in oceanography by Bleck and Boudra (1981) and much developed since then (see Bleck 2002; Hallberg 1995), calculate potential energy in density coordinates. Surprisingly, however, the atmospheric (oceanic) energy equations discussed by Bleck (1978, 1985) in isentropic (isopycnic) coordinates do not show the sought-after expression for the energy equation in density coordinates given below in (18). While Lorenz’s original concept of available potential energy is a global one, Holliday and McIntyre (1981) found an expression for the *local density* of available potential energy for an incompressible stratified fluid, which has become the standard basis of most of the recent developments:

where the reference state is a state of rest given by the density profile *ρ*_{r}(*z*), and the expression gives the local work of the gravity and pressure forces required to move a fluid parcel from the reference state at depth *z** to the actual state at depth *z* along a density preserving path. That is, *ρ*(*x*, *y*, *z*, *t*) = *ρ*_{r}(*z**). For small vertical displacements, a Taylor expansion to first order gives

so that (1) becomes approximately

where *η* = *z*′ − *z****. The integration gives the well-known expression for available potential energy correct for small vertical displacements *z* − *z**:

with the Brunt–Väisälä frequency given by

With the definition of the perturbation density *ρ*′ = *ρ*(*x*, *y*, *z*, *t*) − *ρ*_{r}(*z*), the small vertical displacement and perturbation densities are linked (in the absence of mixing) by

so that (2) can be rewritten as

Relations (2) and (4) are the “textbook” expressions that appear in the context of small-amplitude internal waves (Turner 1973; Lighthill 1978; Gill 1982) and quasigeostrophic theory (Pedlosky 1987). Variants of expression (4) have been used to evaluate eddy potential energy associated with the mesoscale motions in the central waters of the North Atlantic by Dantzler (1977) and the mean World Ocean potential energy by Oort et al. (1989). In either case, the vertical displacements of actual density surfaces may not be small, and Roullet and Klein (2009) have pointed out the important role of higher-order *anharmonic* effects for geostrophic turbulence. Roullet et al. (2014) returned to the original definition in (1) to make an improved estimation of *eddy* potential energy for the World Ocean from data of the Argo float program. Scotti et al. (2006) and Kang and Fringer (2010) showed that approximation (2) is also a poor one for nonlinear internal waves and that the former definition in (1) must be used.

When deriving the expression for potential energy with a Lagrangian density coordinate, we find that the natural quantity to consider is the top-to-bottom vertical integral of the potential energy [i.e., the top-to-bottom vertical integral of (1)]. The momentum equations that are chosen to derive the corresponding energy equation are the traditional primitive equations. The underlying hydrostatic pressure assumption requires a small aspect ratio of vertical scales over horizontal scales, which is well justified in the ocean for the *small slopes* of density surfaces from the mesoscale to the larger scales of the circulation. In water mass formation regions, frontal regions, and high-frequency internal waves, to name a few, they become vertical, and nonhydrostatic effects come into play.

Although likely to be known to some readers, the expression for the potential energy of an entire fluid column is derived in section 2 in terms of the vertical displacements of the density interfaces. Less well known, however, is the energy equation derived in section 3, with a new and direct proof that shows that the potential energy of the shallow water system coincides with that derived from first principles. The determination of the Lorenz reference state is covered in section 4. We apply these results to the simulation of the double-gyre, wind-driven circulation, as calculated with a multilayer shallow water model in section 5, and to the energy determination of the time-mean oceanic circulation (as reconstructed with the help of Argo float mean displacements) in section 6. We observe that the gradient of the time-mean available potential energy becomes a useful predictor for the presence of eddy energy.

## 2. The potential energy of a fluid column with a density coordinate

Montgomery (1937) suggested rather early the use of potential density as a vertical coordinate in studies of ocean circulation, from the perception that the circulation was flowing primarily along potential density surfaces. In many theoretical or modeling works, stratification is discretized as layers of constant density, separated by interfaces assumed to be material surfaces under conservative (adiabatic) physics. Vertical velocities need not be explicitly considered in these multilayer shallow water models, an advantage that has been exploited in many studies: for example, the ventilated thermocline of Luyten et al. (1983), MICOM (Bleck and Boudra 1981), the Hallberg isopycnal model (HIM) and GOLD model (Hallberg 1995; Adcroft and Hallberg 2006), and the Hybrid Coordinate Ocean Model (HYCOM; Bleck 2002). A comparison of the modeling of ocean circulation with depth and isopycnic coordinates can be found in Chassignet et al. (1996).

Consider a layer of density *ρ*_{i} limited above by an interface at *z* = *η*_{i}(**x**, *t*) and below by an interface at *z* = *η*_{i+1}(**x**, *t*), with **x** being the position vector in the horizontal plane. The free surface is at *z* = *η*_{1}(**x**, *t*), the bottom is at *z* = *η*_{B}(**x**), and an arbitrary number of *n* such layers is chosen. The topography is arbitrary. The density of the layers is monotonically increasing since the fluid is stably stratified: *ρ*_{i+1} − *ρ*_{i} > 0, ∀*i* (see Fig. 1). The gravitational potential energy of a fluid particle of mass *m*, density *ρ*, and volume *dz **dA*, with *dA* an infinitesimal horizontal area, is simply *ρ **g **z **dz **dA* with the *z* axis directed upward. The potential energy per unit area of a full water column is then

With the arrangement of layers and interfaces in Fig. 1, this can be written successively as

As the density is now constant in each integral,

Grouping together adjacent terms,

where the contribution of the free surface to potential energy is represented by the first term, the contribution of internal baroclinic displacements by the second, and the contribution of the bottom topography by the last term. Note that an absolute determination of the free surface, hence, an absolute circulation, is required to estimate the first term. The calculation for a continuous (stably stratified) distribution of density proceeds as follows. Equation (5) is first integrated by parts to give

where *η*_{b} (*ρ*_{b}) and *η*_{s} (*ρ*_{s}) are the height (density) of the bottom and free surface, respectively. Making a further change of variables between *z* and *ρ* coordinates, the last term now transforms to an integral with respect to *ρ*:

Equation (8) is readily seen as the continuous version of (6). This last equation is identical with (8) of Oort et al. (1989), apart from the absence of boundary terms. Curiously, Oort et al. (1989) did not use this equation for their estimates of the potential energy of the World Ocean, but proceeded instead with the quasigeostrophic version in (4). One reviewer has pointed out that (6) or very similar APE equations have been used for decades in the diagnostics of several ocean models [Poseidon, Parallel Oregon State University Model (POSUM), HIM, and MOM6, at least]. We now show that (6) is the appropriate form of potential energy for the multilayer shallow water model.

## 3. The energetics of the shallow water model

As far as we know, the first mention of the proper form of potential energy of the multilayer shallow water model is by Ripa (1991) in the context of the derivation of hydrodynamic stability theorems. A detailed derivation of the energy equation of the two-layer version can be found in Cushman-Roisin and Beckers (2011). The objective of the following derivation is to provide a new proof that avoids the Boussinesq approximation and is more direct than Ripa’s derivation.

Inside a layer of constant density in Fig. 1, the horizontal pressure gradient that drives the motion is independent of *z*, and so is the horizontal velocity **u**_{i}. This pressure derivative with respect to *x* (at constant *z*) can be written in density coordinates as

Using hydrostatics , this expression can be rewritten as

and similarly for *y* derivatives, where the new variable *M*, the Montgomery potential, is simply the dynamic pressure:

The inviscid and adiabatic momentum equations for a layer *i* of constant density can now be written in terms of *M* as

In the absence of mixing, the continuity equation is simply

Note that Einstein summation convention does not apply in (9) and (10). The thickness of layer *i* is *h*_{i} = *η*_{i} − *η*_{i+1}, and from now on, the ∇ symbol in (9) and (10) is the 2D horizontal gradient *at constant density* in Cartesian, spherical coordinates, or locally orthogonal coordinates as appropriate. To derive the energy equation, the scalar product of *h*_{i}**u**_{i} with the momentum equation [(9)] is carried out. A few identities are useful:

Similar operations allow us to write

The first term represents the local rate of change of kinetic energy per unit volume; the second is the divergence of the kinetic, potential energy, and pressure flux; and the last term remains to be identified. As it turns out, the vertical sum of this last term is precisely the rate of change of potential energy [(6)]. To show this, the relation between pressure and interface height must be found from hydrostatics (see, e.g., Pedlosky 1987; Cushman-Roisin and Beckers 2011). Leaving aside an atmospheric pressure contribution, the dynamic pressure in layer 1 is simply *p*_{1} = *g**ρ*_{1}*η*_{1}, the Montgomery potential *M*_{1} of layer 1. In the interior, integrating hydrostatics from *z*_{i−1} to *η*_{i} and from *η*_{i} to *z*_{i} allows us to relate the pressure at *z*_{i−1} and *z*_{i}:

Using continuity of pressure at the interface, the sum of these equations becomes

or, in terms of the Montgomery potential *M* = *p* + *ρgz*,

Equation (14) is the form taken by the hydrostatic relation in density coordinates: it relates the jump of the Montgomery potential from one layer to the next to the height of the interface and the density jump across it. To get the total energy equation, we have to sum over the *n* layers. This sum is first expanded as

(the time-independent bottom topography term *η*_{b} drops out). Using (14) and *M*_{1} = *g**ρ*_{1}*η*_{1}, this finally gives

an expression that is nothing but the rate of change of total potential energy density per unit area *P*_{E}:

This expression is identical with (6) (apart from the constant bottom term). Define similarly the total kinetic energy density *K*_{E} (per unit area) as

The layer-wise sum of the energy equation [(12)] gives the full energy equation:

The area integral of (18) over a closed ocean basin with fixed boundaries readily demonstrates energy conservation in the absence of mixing and friction:

Note that the kinetic and potential energy density *K*_{E} and *P*_{E} are expressed here in joules per horizontal area (J m^{−2}), whereas the potential energy density of Holliday and McIntyre (1981), used by MacCready and Giddings (2016), Tailleux (2013), Kang and Fringer (2010), and Lamb (2008), is local in the vertical (J m^{−3}).

## 4. The available potential energy in the multilayer shallow water model

Only the part of potential energy that can be converted to kinetic energy over a given domain is physically significant. On a rotating Earth, the Rossby adjustment problem shows the difficulty of converting potential to kinetic energy: the final geostrophic state still has much potential energy that cannot be further converted [see the discussion in Gill (1982); Cushman-Roisin and Beckers (2011)]. Lorenz (1955) chooses a reference state for potential energy with two properties: 1) it is a state of rest, hence the density interfaces are flat, and 2) the transition from the reference state to the actual state is carried out through adiabatic (conservative) processes. Lorenz’s available potential energy is then the difference between the total potential energy and that of the reference state. The formulation of potential energy in density coordinates is well adapted to the determination of that reference state because Lorenz’s idea leads directly to the conservation of volume of layers of constant density between the actual and the reference states. The *constant* heights *Z* of an arbitrary reference state *at rest* are defined by a certain function *Z*(*ρ*) (to be determined later). Using (16), the potential energy associated with this reference state is therefore

where this reference potential energy is simply a constant independent of horizontal position (the *Z*_{i} and corresponding Montgomery potentials of the reference state are all independent of horizontal position). Therefore, the difference obeys the energy equation [(18)]. Selecting the reference state chosen by Lorenz, only adiabatic, conservative, particle displacements between the reference state and the actual state are allowed. The *vertical displacement ζ*_{i} of a fluid particle in layer *i* at horizontal position **x**_{0}, and vertical position *Z*_{i} in the reference state at *t* = 0, which moves adiabatically (staying on the same side of a density interface) to another horizontal position **x** and vertical position *η*_{i} at time *t*, is *η*_{i} = *Z*_{i} + *ζ*_{i}. Assume first a flat bottom for simplicity, and define the spatial horizontal average of any quantity *s* over the domain by . If a density interface outcrops at the surface, the interface is simply continued at the surface to keep the same number of interfaces at each position. Because the area occupied by any interface is the same, the Lorenz requirement of the conservation of the volume of each layer of constant density is

But since we can always choose the surface geopotential to be such that 〈*η*_{1}〉 = 0, we obtain by recurrence that

The height of any reference state *Z*_{i} is simply the horizontal area average of the actual interface height. Introducing *η*_{i} = *Z*_{i} + *ζ*_{i} in the actual potential energy [(16)] and using 〈*ζ*_{i}〉 = 0 gives

We recognize in the first two terms the potential energy of the reference state [(20)], while the last two define the mean perturbation potential energy (i.e., the mean available potential energy):

Expression (22) shows that the available potential energy can be written entirely in terms of the *variance of the vertical component of the particle displacements from the reference state*. It represents the energy associated with the spatially varying density stratification, the part that may be exchanged with kinetic energy. We further define the unaveraged available *perturbation* potential energy as

This expression (23) is interesting since it allows us to find out the geographic distribution of the local available potential energy in the domain. The vertical displacements in (22) and (23) can be as large as desired, provided the slopes of the isopycnal displacements remain small to stay in the hydrostatic regime. When topography is present and density interfaces intersect the bottom, it is possible to continue to express APE as the variance of vertical displacements from the reference state, provided the missing density interfaces are continued along the topography. By doing so, the reference state can be computed as the spatial average of density interfaces as previously, and APE follows as with (22) and (23). The appendix details these considerations.

## 5. Energetics of the wind-driven shallow water model

As a first step in the investigation of this APE definition and its link with the mean circulation, we have tested it in a traditional double-gyre wind-driven experiment (Bleck and Boudra 1986). We use HYCOM for its parallel computing capabilities, modified with the help of Yves Morel, Eric Chassignet, Alan Wallcraft, and Alexandra Bozec, for implementing the adiabatic shallow water equations and validate the results at low resolution against a purely isopycnal version of MICOM (as modified by Herbette et al. 2003). No mixed layer or diapycnal diffusion is used, and the wind stress is evenly distributed in the upper 50 m. The four-order momentum advection scheme is used with a combination of deformation-dependent and speed-dependent biharmonic viscosity (Winther et al. 2007). The flat-bottom ocean basin is rectangular in Cartesian coordinates, extending from 15° to 55°N (4448 km) and is 60° wide in longitude (6672 km) and 3800 m deep. A constant zonal wind stress forcing is used, varying with latitude only, driving a double gyre with reduced wind in the subpolar region, compared to the subtropical region:

where *θ*_{0} = 38°N, *θ*_{1} = 15°N, *θ*_{2} = 55°N, and the wind stress amplitude *τ*_{0} = 0.06 N m^{−2}.

We use six isopycnal layers with density anomalies of 25.1871, 25.9636, 26.5789, 27.0365, 27.2754, and 27.3370 kg m^{−3} [based on an analytical exponential profile fitted on the globally averaged *World Ocean Atlas*, with initial thicknesses of 200, 250, 350, 600, 1000, and 1400 m; Locarnini et al. (2010)].

A series of experiments was performed from 1/3° to 1/96° horizontal resolution. We show here the results with the 1/48° (2.3 km) resolution, requiring a baroclinic (barotropic) time step of 240 (12) s. The model is integrated from the initial rest state for 100 years, the spin up takes more than 30 years, and statistics are computed over the last 50 years. APE is computed according to (23) and decomposed into temporal mean APE (MAPE), computed from the 50-yr-averaged SSH and layer-thickness fields and eddy APE (EAPE). In the same manner, the kinetic energy *K*_{E} is computed from (17) and decomposed into mean kinetic energy (MKE) and eddy kinetic energy (EKE). In total, 94.5% of the available potential energy is in the temporal mean part, reflecting the large-scale, large-amplitude variations of isopycnals, whereas 90.5% of kinetic energy is in the eddy part (90.5%), reflecting the efficient baroclinic instability at the Rossby radius deformation scale. The mean surface height is shown in Fig. 2f. The western boundary current exits as a zonal meandering jet near 32°N. The jet broadens and moves to 35°N after 20°E. The mean APE shows a two-lobe structure for the maxima in the subtropical and subpolar gyres separated by a valley lying along the meandering jet, which separates the two gyres. The mean APE exceeds the MKE by two orders of magnitude (ratio = 109). Most interestingly, the regional distributions and intensities of EAPE and EKE match closely. Given that the eddies are generated by baroclinic instability, we expect them to have a link with the locally available mean APE. To do so, we have computed the modulus of the gradient of the mean APE in Fig. 2e. The close comparison of the spatial structure of the two forms of eddy energy in Figs. 2c and 2d with the structure of the gradient of mean APE in Fig. 2e shows that the eddies maximize their energy along the strong gradients of mean APE in the western part of the intergyres. Note how these maxima encircle the APE minimum seen between the two gyres on Fig. 2a. This result in the simplest wind-driven double-gyre setup confirms the interest of such APE calculation following the energy conservation equations. The next step is to find out if this result holds in the realistic global ocean.

## 6. Mean available potential energy, mean and eddy kinetic energy in oceanic basins

In the following, we estimate the MKE and available potential energy of the observed time-mean circulation using a reconstruction of the circulation from the *World Ocean Atlas* (Antonov et al. 2010; Locarnini et al. 2010) and the Andro database of Argo floats displacements at 1000 dbar (Ollitrault and Rannou 2013). Interesting connections of the mean APE and surface EKE as obtained from satellite altimetry are found.

### a. Mean kinetic energy

We first estimate mean kinetic energy using traditional pressure coordinates. Ollitrault and Colin de Verdière (2014) have calculated an absolute geopotential at 1000 dbar from the time-mean motions of Argo floats at that level. Using in situ density calculated from the *World Ocean Atlas* (Antonov et al. 2010; Locarnini et al. 2010) at the 1° scale, we use hydrostatics to estimate an absolute geopotential from the top to the bottom of the ocean. Geostrophy then provides the horizontal velocity and horizontal kinetic energy. Integrating vertically provides the local density of mean kinetic energy (Fig. 3). MKE is large in the western boundary currents (North Atlantic, North Pacific, and South Indian), in the Agulhas retroflection from the tip of South Africa to the Kerguelen Plateau, and in the Southern Ocean. In the Antarctic Circumpolar Current (ACC), mean kinetic energy is organized as elongated patches oriented in the northwest–southeast direction. The unique Zapiola Gyre stands out as a solitary feature at 40°S, 40°W. The zonal tongues of high kinetic energy in the tropical Pacific are conspicuous. Away from these energetic regions, the mean kinetic energy drops by two to three orders of magnitude over vast interior regions.

### b. Mean available potential energy

The fundamental issue is the choice of the density interfaces, which are assumed to be material surfaces under conservative physics. Potential density was used by Oort et al. (1989), but its nonmonotonic character in some regions is not suitable for a vertical coordinate. Several developments have appeared since then: the orthobaric density of de Szoeke et al. (2000), the neutral density introduced by Jackett and McDougall (1997), and the generalized patched potential density and the thermodynamic neutral density of Tailleux (2016). These different quantities attempt to guarantee the adiabatic character [i.e., material conservation and the neutral character that is absence of buoyancy force along particle paths; see Tailleux (2016) for an examination of these properties]. The neutral density *γ*_{n} computed by Jackett and McDougall (1997) is chosen here because it is adapted to the present climatological data and is the variable with the most neutral property. Its nonconservative character is discussed by McDougall and Jackett (2005), and Tailleux (2016) shows that neutral density can be approximated very closely by the conservative thermodynamic neutral density. The temperature and salinity data of the *World Ocean Atlas 2009* (Antonov et al. 2010; Locarnini et al. 2010) and the Jackett and McDougall (1997) algorithm allow us to compute neutral density at each location of the 1° × 1° grid. A simple linear interpolation is then used to find the interface depth *η*_{i} of each of the above *γ*_{ni} values when they exist. When they are not defined, the density interfaces are continued at the surface or at the bottom (see appendix). We have chosen 17 neutral density values *γ*_{n} as vertical coordinates from 24.5 to 27.95, whose distribution remains monotonic at all points of the domain (see Table 1). We have not included any larger*γ*_{n} values to resolve the varieties of bottom and deep waters present in semienclosed basins. The largest density value *γ*_{n} = 27.95 lies at a mean depth close to 2000 m, with maximal excursions reaching 5000 m. The reference state for the available potential energy is found by assuming global volume conservation of the layers between adjacent density interfaces. When density interfaces disappear, they are continued along the topography or along the surface; the depths *Z*_{i} of the reference state shown in Fig. 4 become simply the area average of the height of the density interfaces as in the flat bottom case. To compute the potential energy, the density jumps across the interfaces that appear in the hydrostatic equation in (14) must also be defined. In observations or models, the pressure is computed from the in situ density and the thermal wind from the horizontal gradient of in situ density. However, as the differences between the horizontal gradients of neutral density and of in situ density are small and consistent with our definition of the interfaces and the reference state, we have made the choice here to use neutral density throughout. To conserve the overall density content between the reference state and the actual state, the density of any layer *i* has been computed by the volume average of all neutral density values that happen to fall between the upper *η*_{i} and lower interface *η*_{i−1} of that layer. These layer densities are given in Table 1. The absolute geopotential computed previously allows us to find the free surface height *η*_{1} (referred to *Z*_{1} = 0) for the first term in (23). The computation of available potential energy from (23) is then immediate and shown in Fig. 5. The tropics are regions of low APE that broaden to the east. Two major contributions can be identified in the figure: the contribution of the subtropical gyres, associated with the depressions of midthermocline isopycnals (deep light fluid), and the deep water contribution, associated with the largest density values found near the surface along the Antarctic continent and the Labrador Sea (shallow heavy fluid). The latter contribution dominates the APE density. The vertical and latitudinal structures of mean kinetic energy and APE are shown in Figs. 6 and 7,respectively. The mean kinetic energy is strongly trapped at the surface and decreases monotonically below 1000 m. The World Ocean, eddy-resolving model of von Storch et al. (2012) reproduces this trapping of MKE above 1000 m but also shows a minimum near 3500 m, which is not found in the observations. The ubiquitous surface trapping of mean kinetic energy can be seen as a confirmation of the assumption of Wunsch (1996) to minimize kinetic energy at the bottom in inverse models of the circulation. It is only in the Antarctic Circumpolar Current near 47°S that a significant amount of kinetic energy is observed to reach the bottom. The vertical structure of terms of the sum [(23)] for mean APE is shown in Fig. 6 (area average) and Fig. 7 (zonal average). In contrast to kinetic energy, the APE increases with increasing *γ*_{n} values. The important contribution of the bowls of the subtropical gyres is found between *γ*_{n} values of 26 and 27.5 and from 20° to 40° latitudes. Large APE values of the deep water are associated with the deepest interfaces rising by more than 2000 m to outcrop in the vicinity of the Antarctic continent and the Labrador Sea. Note the very small APE values for low *γ*_{n}, but there is a contribution spanning the tropics around *γ*_{n} = 24.5, which comes from the Pacific.

Using quasigeostrophic theory, Gill et al. (1974) showed that the ratio KE/APE should scale as the Burger number, given by the ratio Rd^{2}/*L*^{2}, with Rd being the first Rossby radius of deformation and *L* the scale of the circulation. However, the geographic dependence of the ratio KE/APE as the square the internal Rossby radius of deformation is just not observed, presumably because the oceanic circulation selects a different scale *L* at different latitudes. When APE and KE are integrated over the World Ocean [from 64.5°S to 64.5°N and excluding all marginal seas, to compare with Oort et al. (1989)], a total mean available potential energy of 2.03 × 10^{20} J is obtained. The deep water contribution explains 46% of this total (isopycnals between 27.7 and 27.95), while 42% originates from the subtropical gyres (isopycnals between 26.5 and 27.6). The contribution of the free surface is only 0.2% of the total. The total mean kinetic energy amounts to 3.41 × 10^{17} J, nearly three orders of magnitude smaller with a global ratio APE/KE near 600. Gill et al. (1974) used a simple 1½ model of a midlatitude thermocline forced by Ekman pumping to find a ratio APE/KE = 10^{3}, quite a remarkable prediction in retrospect, given the simplicity of their model. The comparison in Table 2 shows that the present results are close to the quasigeostrophic estimates of Oort et al. (1989) and Peixoto and Oort (1992). Other estimates from observations using different methods (Huang 2005) or data assimilation models (Zemskova et al. 2015) are nearly an order of magnitude larger.

### c. The connection with eddy kinetic energy

Because results from realistic models indicate that eddy energy originates from the APE reservoir (e.g., von Storch et al. 2012), we asked whether a visible, local connection between the two could be identified. The idealized HYCOM results of the previous section give the clue that eddy kinetic energy develops along the gradients of APE, and the objective is to test this prediction in observations. The observed oceanic geostrophic turbulence is now known globally from surface drifters (Lumpkin and Johnson 2013), from surface EKE from satellite altimetry (Dibarboure et al. 2011), from eddy kinetic energy at 1000 m obtained from Argo float displacements (Ollitrault and Colin de Verdière 2014), and from eddy potential energy also from Argo float profiling *T*/*S* capabilities (Roullet et al. 2014). As these different products present a high degree of spatial coherence at midlatitudes, we use the eddy kinetic energy that has been recalculated from geostrophy using the altimetric data from Dibarboure et al. (2011). First of all, the spatial coherence of eddy kinetic energy and mean kinetic energy confirm the early indications of Stammer (1997). Because the exchange between eddy energy and APE is not expected to have a global character (an eddy in the Gulf Stream hardly depends on APE in the Pacific), we focus on regional basins and compute APE with a local Lorenz state calculated for that basin. If eddies get their energy from APE, it is not expected to observe high EKE at locations of high APE. It is more likely to observe high EKE near an APE minimum or along the flanks of the APE distribution, and we propose to compare EKE with the *gradient* of APE. This gradient (more precisely, its modulus) was calculated with centered differencing of APE summed over density intervals associated with the potential energy stored in the bowls of the main thermocline. We compare such gradient of APE with surface eddy kinetic energy for the Gulf Stream, the Kuroshio, the Agulhas retroflection, and the confluence region.

Starting with the Gulf Stream (Fig. 8), two areas of maximum gradient of APE are found from 70° to 45°W between 35° and 40°N, separated by a minimum value of APE. This minimum is found where the reference state intersects the strongly sloping density interface of the Gulf Stream, a situation similar to the model result in Fig. 2. Note simply the corresponding maxima of EKE found somewhat downstream around 60°W. This geography is consistent with the idea that the gradients of APE locate the source of the eddies that grow and are advected downstream by the Gulf Stream. This correspondence of the two fields continues far in the interior along the Azores Current and within the North Atlantic drift. Note also how small gradients of APE also coincide with low eddy kinetic energy in the eastern part of the basin. A similar situation is found in the Kuroshio region (Fig. 9). The agreement is quite good in the western subtropical gyre, but the correspondence is less satisfying in the northeast: there is eddy energy along a narrow band along 50°N, but it is low otherwise. Yet, a significant gradient of APE is found in a broader region, whose origin comes from the APE stored in the outcrops of density interfaces. In the Agulhas retroflection (Fig. 10), the eddy kinetic energy peaks around 20°–30°E and is associated with the presence of the minimum of APE surrounded by the two maxima of the gradient of APE to the north and south. The eddy energy that leaks into the South Indian Basin follows nicely the contours of the gradient of APE. As in the Pacific, there is little eddy energy linked to the southern maximum of the gradient (caused by the outcrops). The last example is the confluence region, one of the most turbulent regions of the World Ocean (Fig. 11). The circulation offshore is dominated by the anticyclonic Zapiola Gyre. The contours of the maximum gradient of APE and of maximum eddy energy are found to encircle the Zapiola Gyre. Eddy energy is western intensified and advected by the anticyclonic mean circulation of the Zapiola Gyre.

These four examples support the idea that the eddy kinetic energy is found in regions of large gradient of APE. Near the western boundaries, this region encircles a thin region of minimum APE. Although the APE is the energy reservoir for midlatitude eddies, the eddies are associated with the regions of fastest growth of unstable waves, and a necessary condition for instability is that the potential vorticity (PV) gradient along isopycnal surfaces changes sign (Charney and Stern 1962). One therefore expects this necessary condition to be fulfilled in the source region of the eddies. Figure 12 shows the potential vorticity for the North Atlantic on the isopycnal *γ*_{n} = 27.1, an isopycnal surface that is at 400-m depth at the eastern boundary and reaches 750 m south of the Gulf Stream (not shown). The relative vorticity term is not included because it makes a small contribution at this resolution. The potential vorticity *f*/*h* increases from tropical latitudes, but there is a maximum in the recirculation region, and a tongue of low-PV water is found just north of it, along the axis of the Gulf Stream. The meridional PV gradient vanishes there, and the existence of this low-PV tongue allows the necessary condition for instability to be satisfied. The curves of vanishing PV gradient and areas of large APE gradient both indicate the same high eddy energy regions. A similar configuration is found at a somewhat shallower level, *γ*_{n} = 26.5 in the Kuroshio. The situation in the Agulhas Basin is different. The instability there appears to be associated with an anomaly of low PV just south of Cape of Good Hope, which is connected with the Antarctic Circumpolar Current. In the confluence basin, there is also a low-PV tongue in the region of the previously mentioned high gradient of APE.

The presence of eddy kinetic energy in the ocean has long been rationalized through large growth rates computed from local baroclinic instability calculations. Stammer (1997), Smith (2007), and Tulloch et al. (2011) point out the geographic similarities of high eddy energy, fast growth rate of baroclinic instabilities, and Eady growth rate . These fast growth rates are therefore associated with strong vertical shear of mean currents, hence high KE. The connection of these fast growth rates with the gradient of APE can be seen as follows. Take the modulus of the gradient of (23):

The above expression can now be transformed under the quasigeostrophic approximation. Since the large-scale currents are geostrophic to a good approximation, we have

Substituting the velocities for the slopes of the density surfaces gives

an expression where the Boussinesq approximation has been made as well (*ρ*_{i} ≈ *ρ*_{0}). The vertical shear signature of the baroclinic terms in the above sum coincides with the dependence of the Eady growth rate.

From these energy considerations, we add the following energy viewpoint to the regimes identified in the early analysis of satellite observations by Stammer (1997):

The tropical regime shows low APE and small gradients thereof, hence little potential energy available for conversion.

The central and eastern energy deserts of the low-latitude subtropics with low EKE, KE, and small gradients of APE. This is the region where linear planetary waves have been found (Chelton and Schlax 1996).

The western boundary current regimes where KE and EKE are high. The gradient of APE diagnostic that we propose indicates the regions of exchange between APE and EKE in the Gulf Stream, the North Pacific, the Agulhas retroflection, and the confluence region (and the northern flanks of the ACC; not shown).

## 7. Summary

The available potential energy in density coordinates has several useful properties, summarized here:

The expression for the vertically integrated APE has the character of a quadratic form (i.e., a weighted vertical sum of the squares of the vertical displacements of density surfaces), which is valid for arbitrary vertical displacements (at the difference of the quasigeostrophic expression).

This expression is the appropriate form for energy conservation in the multilayer shallow water model.

The application of this formulation in the context of data originating from

*z*coordinate models or in situ observations is straightforward, requiring simply an interpolation procedure to determine the density interfaces. The conservation of volume between any two interfaces implies that the Lorenz reference state reduces to the area average of the interface heights.

The local formulation of Holliday and McIntyre (1981) is also valid for arbitrary vertical displacements but does not share the quadratic property. Indeed, Roullet and Klein (2009) draw attention to anharmonic terms (beyond quadratic) that appear for large vertical excursions when the Holliday and McIntyre (1981) local energy density is used. The present formulation allows, in turn, insights into the stability properties of the mean oceanic circulation: the high-resolution simulations show that the gradient of mean APE predicts well the gross region of large eddy kinetic energy. We have also applied our method to find the APE of the World Ocean, using the neutral density as an approximation to the material conservation property, which is required to find the Lorenz state. Further improvements in our method are linked closely to that issue. It is planned in the near future to carry out a comparison of our method with Tailleux (2013) and Saenz et al. (2015) using a common dataset. We have compared APE with the mean kinetic energy and found that the former is nearly three orders of magnitude larger than mean kinetic energy, quantifying the small overall Burger number of the time-mean oceanic circulation. The encouraging comparisons of eddy kinetic energy obtained from satellite altimetry with the regions of large gradient of mean APE encircling an APE minimum in western boundary current regions confirms the HYCOM results. The diagnosis appears also consistent with the local necessary condition of baroclinic instability (a reversal of the sign of the potential vorticity gradient along isopycnal surfaces) and the *local* baroclinic instability calculations (Smith 2007; Tulloch et al. 2011). Of course, each region will require a far more careful examination than has been provided here.

The determination of the reservoirs of available potential energy and kinetic energy in the real ocean is going to improve regularly with the growing extent of the Argo float program, but it is impossible, for the time being at least, to obtain the fluxes between mean APE and eddy energy since they require the knowledge of vertical velocities at the eddy scale. Given this difficulty, the search for correspondence between eddy energy and gradients of mean APE in both eddy-resolving model simulations and oceanic observations is worth exploring further.

## Acknowledgments

We thank E. Chassignet, A. Hochet, Y. Morel, F. Roquet, J. Sirven, and R. Tailleux for comments on an early draft of this manuscript. The manuscript was much improved by the comments and questions of the reviewers, P. MacCready and an anonymous one. They also pointed out an inaccurate treatment of the topography in the World Ocean APE calculation. The authors acknowledge the Pôle de Calcul et de Données Marines (PCDM at Ifremer, Brest) for providing DATARMOR computational resources for the HYCOM simulations.

### APPENDIX

#### Calculation of APE with Arbitrary Topography and Density Interfaces

The calculation of the depths of the Lorenz state and APE is now described for the general case of varying topography. We illustrate the problem with a three-layer model described in Fig. A1 that captures all the difficulties that have to be taken care of. The lower interface *η*_{2} of layer 1 outcrops at point A. When it is not defined, this interface is now carried on at the surface (on the left of point A) with *η*_{2} = *Z*_{1} = 0. Similarly, when the lower interface *η*_{3} of layer 2 incrops, such as at point B, the interface is carried on at the depth of the bottom; hence, *η*_{3} = −*H*. This continuation allows us to compute the volume of each layer over the *same domain of integration* irrespective of the existence of a given interface:

Since the horizontal area *A* is now the same for all layers, the thicknesses are just *H*_{i} = *V*_{i}/*A*, and starting from *Z*_{1} = 0 at the surface, the height *Z*_{i} is found from *Z*_{i} = *Z*_{i−1} − *H*_{i−1}. With the continuation procedure, the height *Z*_{i} becomes, quite simply, the spatial average 〈*η*_{i}〉, exactly as for the previous expression for flat bottom and no outcropping–incropping of density interfaces.

The last thing is to find out how to compute available potential energy. In formulations (22) and (23), the last topographic term present in (6) is absent since the bottom density is constant. But this is not the case anymore. It is possible, however, to keep the previous formulation in the general case, as the extreme frontal situation of Fig. A2 shows. Two water masses of equal volume, equal area Σ, and different density now lie side by side in configuration A. The interface between the two layers is continued at the bottom *η*_{2} = −*H* and at the surface *η*_{2} = 0. The two volumes are then readjusted to the Lorenz reference state in B, where *Z*_{2} = −*H*/2. With this definition of *η*_{2}, expression (16) shows that the potential energy for each sides of the domain can be written as

And the total APE per unit area recovers the well-known result

But the APE can also be computed more symmetrically by introducing the height anomaly *ζ*_{2} = *η*_{2} − *Z*_{2}, which is −*H*/2 on the left and +*H*/2 on the right:

with the same result as before for the total APE. This last result shows that expressions (22) and (23) can still be used in the general case, provided the APE calculation includes also the surface and bottom contributions of the stacked density interfaces. The topographic case can be handled in the same way. Including a step of depth *H*_{1} (<*H*) for the box of density *ρ*_{1}, the answer is provided by the same APE formula with *H* replaced by *H*_{1}. The APE of the denser water in the deeper box 2 vanishes since its volume remains at the same position in the reference state.

## REFERENCES

*Salinity*. Vol. 2,

*World Ocean Atlas 2009*, NOAA Atlas NESDIS 69, 184 pp.

*Introduction to Geophysical Fluid Dynamics: Physical and Numerical Aspects*. 2nd ed. International Geophysics Series, Vol. 101, Academic Press, 828 pp.

*Atmosphere-Ocean Dynamics*. Academic Press, 662 pp.

*Waves in Fluids*. Cambridge University Press, 504 pp.

*Temperature*. Vol. 1,

*World Ocean Atlas 2009*, NOAA Atlas NESDIS 68, 184 pp.

*Geophysical Fluid Dynamics*. Springer-Verlag, 710 pp.

*Physics of Climate*. Springer, 520 pp.

*Buoyancy Effects in Fluids*. Cambridge University Press, 368 pp.

*The Ocean Circulation Inverse Problem*. Cambridge University Press, 442 pp.

## Footnotes

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).