Eddy Cancellation of the Ekman Cell in Subtropical Gyres

The presence of large-scale Ekman pumping associated with the climatological wind stress curl is the textbook explanation for low biological activity in the subtropical gyres. Using an idealized, eddy-resolving model, it is shown that Eulerian-mean Ekman pumping may be opposed by an eddy-driven circulation, analogous to the way in which the atmospheric Ferrel cell and the Southern Ocean Deacon cell are opposed by eddy-driven circulations. Lagrangian particle tracking, potential vorticity ﬂuxes, and depth–density streamfunctions are used to show that, in the model, the rectiﬁed effect of eddies acts to largely cancel the Eulerian-mean Ekman downwelling. To distinguish this effect from eddy compensation, it is proposed that the suppression of Eulerian-mean downwelling by eddies be called ‘‘eddy cancellation.’’


Introduction
There exist a number of circulations in the ocean and atmosphere that are prominent in Eulerian-mean velocity fields but are inconsistent with observed tracer distributions, for example, the Deacon cell in the Southern Ocean (Döös and Webb 1994). When the effect of eddies is included in the averaging operator, these circulations are dramatically altered or disappear entirely.
Large-scale sinking due to Eulerian-mean Ekman pumping is believed to dominate the vertical motion in subtropical gyres (see, e.g., Oschlies 2002;Williams and Follows 2011) and to form part of an intergyre overturning circulation known as the Ekman cell. A schematic of the Ekman cell is shown in Fig. 1. The diabatic circulation of the Ekman cell is difficult to reconcile with the adiabatic nature of the ocean interior. Cancellation of large-scale downwelling in subtropical gyres by an eddy-induced circulation would remove the substantial diapycnal transport across the subtropical thermocline associated with the Ekman cell. As described by Griffies et al. (2015), if the ocean interior is adiabatic and temperature surfaces are approximately aligned with density surfaces, then we expect the net vertical transfer of heat to be small and the downward transport of heat by the mean flow to be opposed by vertical eddy transports.
Downwelling in the subtropical gyres is also believed to remove nutrients from the oligotrophic surface waters (Williams and Follows 2011). Cancellation of the Ekman cell by eddies may help settle the long-standing discrepancy between new production and estimated nutrient supply in the gyre-scale nutrient budget (see, e.g., McGillicuddy 2016; Follows 1998, 2011). Because such a large proportion of the ocean surface is covered by subtropical gyres, these regions contribute approximately half of all marine export production (Jenkins and Doney 2003;Williams and Follows 2011). This large contribution means that improving our understanding of vertical transports in subtropical gyres has implications for the global nutrient and carbon budgets. The presence of an eddy-induced circulation balancing time-mean Ekman downwelling also has implications for mode-water formation and the maintenance and structure of the thermocline (see, e.g., Henning and Vallis 2004;Marshall et al. 2002;Polton and Marshall 2003;Radko and Marshall 2004b;Xu et al. 2016).
This paper is organized as follows: Section 2 describes the atmospheric Ferrel cell and the Deacon cell in the Southern Ocean, their sensitivity to the choice of averaging technique, and some theoretical background regarding different methods for averaging flows. Section 3 contains the hypothesis that we seek to test in this paper. In section 4, we present the model used for our idealized simulations. In section 5, we present the results of our Lagrangian particle tracking experiments. In section 6, we use vertical fluxes of potential vorticity to examine the balance between vertical potential vorticity fluxes from Eulerian-mean advection and eddies. Section 7 quantifies the extent to which eddies oppose the Eulerian-mean circulation through the use of depth-density overturning streamfunctions. In section 8, we present our concluding remarks and some possibilities for future work that could expand upon our current results.

Averaging flows and vanishing circulations
To develop an understanding of the chaotic and turbulent patterns of circulation commonly observed in the ocean and atmosphere, we need to take an average of the flow. However, the way in which we take this average can have a dramatic impact on the structure, strength, and even direction of the diagnosed circulation (see, e.g., Plumb and Mahlman 1987). In turbulent flows, such as those analyzed in geophysical fluid dynamics, the transport of mass and tracers due to fast, small-scale, timevarying features of the flow field can be substantial. If our averaging operator does not account for these transports, then there may be substantial differences between the transport of the averaged flow and the averaged transport of the flow. While taking an average cannot create a flow where none exists, it can provide a misleading description of the flow if considered in isolation.
The rectified effect of eddies is a combination of downgradient diffusion along isopycnals and advective transport (Gent et al. 1995;Plumb 1979). In certain geometries, this advective transport can lead to upgradient transport (Lee et al. 1997) and can oppose circulations in the Eulerian-mean fields. There are multiple studies showing that correctly representing eddies is a necessary condition for accurately modeling the distribution of tracers. For example, Follows and Marshall (1996) modeled the distribution of a transient radioisotope of carbon and found that models without mesoscale eddies were unable to reproduce the observed distribution.

a. Averaging flows
We restrict our attention to two averaging techniques because they are conceptually simple and encompass the range of behaviors that we wish to discuss. For a more rigorous and mathematically demanding approach to averaging the equations of motion, the reader is referred to Young (2012) and Maddison and Marshall (2013).  . The thick black lines show the circulation of the Ekman cell with fluid sinking in the subtropics and returning to the surface in the subpolar and tropical latitudes. The white arrows show the direction of Ekman transport and pumping, while the gray arrows show the eddy transports that oppose the Eulerian-mean Ekman transport.

1) EULERIAN MEAN
Taking the Eulerian mean of a field is straightforward: a time averaging operation is performed independently at each spatial location. A variable X is separated into a timemean component, usually denoted X, and time-varying component, denoted X 0 . The effect of eddies and other variability is taken into account by calculating the cross correlations between terms. For example, if we take the inviscid Boussinesq zonal momentum equation and apply a time averaging operator, the equation becomes in which u is the zonal velocity, u is the three-dimensional velocity field, f is the Coriolis parameter, y is the meridional velocity, r 0 is the reference density, p is the pressure, x and t are the zonal and temporal coordinates, and the final term represents the effect of the variability on the time-mean field, known as the Reynolds stress. We retain the time derivative of the time-mean zonal velocity, since this term describes changes over time scales comparable with, or longer than, the time scale of the time-mean operator. The Eulerian-mean formalism splits the tendency terms into contributions from time-mean and timevarying components. It is this artificial separation that lead Dunkerton (1980, p. 392) to state that the Eulerian mean ''can be and has been misleading.'' Within the Eulerian-mean framework it is possible to identify circulations that are incompatible with observed tracer distributions, such as the Deacon cell in the Southern Ocean (Danabasoglu et al. 1994;Döös and Webb 1994;. The location of inconsistencies between Eulerian-mean circulations and tracer distributions highlights areas where the time-varying contribution is nonnegligible. Another example of a situation in which the Eulerianmean velocity field is misleading may be found by considering the propagation of planar waves across the surface of a fluid otherwise at rest. The fluid velocity at the wave crests is in the direction of wave propagation, and the particles are lifted and moved forward as a crest passes their position. As a trough passes their location, the particles are lowered and advected in the opposite direction. The magnitude of the wave velocity decays with depth, which leads to the particles experiencing a net displacement after each wave period, despite the fact that below the level of the wave troughs the Eulerianmean velocity is identically zero. The difference between the Eulerian-mean velocity and the velocity following a particle is known as the Stokes' drift velocity. For further discussion of Stokes' drift the reader is referred to chapter 3 of Williams and Follows (2011).

2) GENERALIZED LAGRANGIAN MEAN
The mathematics of the generalized Lagrangian-mean (GLM) framework is formidable, but in this case we require only an intuitive understanding of its implications. GLM was developed as a framework to analyze wave-mean flow interactions in the stratosphere (Andrews and McIntyre 1978a,b). The GLM velocity is not defined following a single fluid parcel but rather following the center of mass of a collection of fluid parcels that span the region being averaged (Andrews and McIntyre 1978a). For example, to calculate a zonal average, a tube of particles would be initialized along a line of latitude. However, after some finite period of time, the particle distribution can be sufficiently altered by the fluid motion that interpretation of the velocity and location of the centroid of the distribution becomes difficult. Figure 2 shows this problem schematically. In the initial circular configuration, the location and momentum of the centroid may be easily associated with the average motion of the surrounding fluid. However, once advection has sufficiently distorted that distribution, shown schematically in the lower portion of Fig. 2, it is no longer clear that the GLM over the parcel is a meaningful quantity. The issue of reinitialization within GLM theory is briefly discussed by Dunkerton (1980), who recommends a case-by-case approach.

1) FERREL CELL
The Ferrel cell is an atmospheric circulation with air rising in the subpolar latitudes and sinking in the subtropical latitudes. Because this circulation involves colder air rising and warmer air sinking, it cannot be sustained by density differences due to temperature. Circulations such as this are known as thermally indirect and require the input of energy to maintain (Vallis 2006). The presence of a Ferrel cell in atmospheric circulations is dependent on the averaging process used to analyze the atmosphere (Dunkerton 1978). Plumb and Mahlman (1987) provide an excellent illustration of how sensitive the existence of the Ferrel cell is to averaging. Plumb and Mahlman present time and zonal means of the velocity fields from a general circulation model in which there are three circulation cells in each hemisphere. However, when Plumb and Mahlman compute the effective transport circulation, which is similar to the GLM but allows for spatially inhomogeneous diffusion, the circulation pattern shows only two cells in each hemisphere, and the Ferrel cell is no longer present. Furthermore, when comparing the two transport fields, they give opposite results for the direction of meridional transport at the surface at a number of latitudes.

2) DEACON CELL
The Deacon cell is a meridional circulation in the Southern Ocean characterized by upwelling of deepwater south of the Antarctic Circumpolar Current (ACC), northward movement of this fluid in the surface Ekman layer, and downwelling north of the ACC (see, e.g., Speer et al. 2000). A Deacon cell can be clearly seen in the timeand zonal-mean circulations of the coarse-resolution ocean model used by Danabasoglu et al. (1994). When Danabasoglu et al. include an eddy-induced bolus velocity in their analysis, there is a near-complete cancellation of the Deacon cell. While most models do not reproduce such near-complete cancellation, they do show a substantially smaller residual overturning circulation that is directed along isopycnals and consistent with the observed water masses Marshall and Radko 2003). The existence and strength of the Deacon cell, which are dependent not only on the flow field but also the choice of averaging operator, are relevant to our understanding of subduction, the formation of deep water, and to tracer advection (see, e.g., Abernathey et al. 2011;Dufour et al. 2012;Marshall and Radko 2003).
The similarities between the Ferrel cell in the atmosphere and the Deacon cell in the Southern Ocean have been explored previously. For example, Karoly et al. (1997) showed that both the Ferrel and Deacon cells are the result of correlations between zonal variations in density and meridional flow, not zonal mean buoyancy forcing. This means that if the averaging is performed using density as the vertical coordinate, neither circulation cell is apparent. Several other studies have also found that the use of density as the vertical coordinate for their time-and zonal-mean averaging operations resulted in the disappearance of the Deacon cell (Döös and Webb 1994;Hirst et al. 1996).

3) EKMAN CELL
The Ekman cell is a diabatic circulation hypothesized to link Ekman-driven downwelling in subtropical gyres with Ekman upwelling in the tropics and subpolar regions. Figure 1, which has been adapted from , shows a schematic of the Ekman cell. The concept of Ekman cells is used extensively in the biogeochemistry literature. For example, Najjar and Keeling (2000) describe a fluid parcel traversing the Ekman cell and use the changes in nutrient concentration and light availability to explain observed productivity and oxygen fluxes. Even when the Ekman cell is not invoked, the surface Ekman flows and the resulting upwelling and downwelling are often used when analyzing the nutrient budget of the subtropical regions (e.g., Williams and Follows 1998).
The geometry of subtropical gyres is more complex than the near-zonal flows in the atmosphere and Southern Ocean. The presence of nonlinear, inertial recirculations and variation in gyre shape with depth makes it more difficult to directly diagnose and analyze the Ekman cell. Polton and Marshall (2007) use horizontal integrals of time-mean vertical potential vorticity fluxes to explore the Ekman cell in a 1 /48 horizontal resolution ocean general circulation model. They find that the time-mean advective flux is partially opposed by an eddy-induced flux and suggest that a higher-resolution model with a better resolved mesoscale eddy field may lead to a more complete cancellation.

Hypothesis
As described above, there are multiple thermally indirect circulations in Eulerian-mean fields that are largely opposed by eddy-driven circulations. Given that the Ekman cell shares many characteristics with these circulations, we seek to test the following hypothesis: there exists an eddy-driven circulation that opposes the Eulerian-mean Ekman cell.

Model
Our model is designed to explore the dynamical effects of resolved mesoscale variability. As such, we focus on resolution, not realism, and have simplified both the geometry and the thermodynamics of our model. The model domain covers 658 of latitude and 308 of longitude. The basin has a continental shelf in the west that is 58 wide and slopes gently downward until reaching 200-m depth, at which point it drops steeply to 4000-m depth. The rest of the domain is a uniform depth of 4000 m. Our model employs a single variable linear equation of state that depends only on temperature. A sinusoidal wind forcing is used to drive the horizontal circulation. We drive a meridional overturning circulation with relaxation of the surface temperature field to a linear function of latitude and relaxation to an exponential temperature profile in a sponge region occupying the southernmost 58 of the domain. The sponge region does not alter any other model variables. The model domain and momentum forcing are shown in Fig. 3.
We use the MITgcm ) in a hydrostatic configuration to numerically solve the equations of motion on a spherical grid with a horizontal resolution of 1 /88 and 40 unevenly spaced vertical levels. The vertical grid spacing is 10 m at the surface and increases to 360 m at the bottom. Midlatitude mesoscale variability is well resolved by our model, as evidenced by the strong, inertial recirculation regions surrounding the zonal extension of the western boundary currents shown in Fig. 3. The model is spun up using a coarser grid, interpolated to 1 /88, and integrated for 60 model years until statistically steady. The simulation is then integrated for a further 10 model years, and the analysis presented here is based on the final 10 years. Momentum dissipation is provided by a combination of Laplacian and biharmonic Smagorinsky viscosity schemes (Griffies and Hallberg 2000;Smagorinsky 1963). The combination of Laplacian and biharmonic viscosity allows lower values for each to be used, resulting in a less viscous simulation (Chassignet and Marshall 2008). FIG. 3. Schematic of the model used in this paper. The surface shows the Eulerian-mean temperature field and Bernoulli potential contours. The eastern face shows the zonally integrated depth-latitude overturning circulation calculated from the Eulerian-mean fields. The southern face shows the continental shelf and rise. The zonal wind forcing is shown to the west of the model domain. Note that the model uses a spherical polar grid but is shown here as a cuboid for simplicity.
The simulation employs the Prather advection scheme (Prather 1986), and thus we expect very little spurious diapycnal mixing (Hill et al. 2012). We include some horizontal and vertical diffusion, with diffusivities of 100 and 10 25 m 2 s 21 , respectively, to reduce grid-scale noise in the temperature field. The model does not use a convective parameterization, such as KPP (Large et al. 1994) or otherwise, because the nonlocal buoyancy fluxes and elevated diffusivities would be problematic for our potential vorticity flux budgets and particle tracking experiments as well as adding to the computational expense. Since our interest lies below the mixed layer in the region of the subtropical mode water and thermocline, we believe this to be an acceptable idealization.

Lagrangian particle tracking
To compute the combined effect of eddies and the mean flow, we use the GLM framework. However, because the shape of wind-driven gyres varies substantially with depth, we use Lagrangian particle tracking techniques to compute an approximate GLM velocity numerically. As described in section 2, the GLM velocity is the mean velocity following a collection of parcels that span the region being averaged. We estimate the GLM velocity in the subtropical mode water of our model by seeding a cloud of particles there. The subtropical mode water in our model is the mass of weakly stratified fluid above the thermocline in the subtropical gyre. We deal with the issue of reinitialization, mentioned in section 2 and discussed by Dunkerton (1980), pragmatically. Our Lagrangian particles are observed to begin leaving the region of interest after 3 months of model integration. Thus, we track particles for 3 months before calculating the center of mass of the particles' distribution and reseeding new particles within an ellipsoid around this point. Only the location of the ellipsoid changes; its size remains constant throughout the experiments.
We use a fourth-order Runge-Kutta scheme to numerically solve the three coupled ordinary differential equations for each particle. Velocity fields are linearly interpolated in space and time to each particle's location. The code we use is part of a Python module for analyzing MITgcm simulations that can be found online (at http://doddridge.me/publications/dmh2016). To test for sensitivity to the time step chosen in our algorithm, we computed a test set of pathways using a time step of 1 min. Using any time step up to 1 h gives results that are indistinguishable for 10 yr of model time; hence, we use a time step of 1 h.
Similarly to Pennel and Kamenkovich (2014), we perform two different particle tracking experiments and artificially alter the velocity field used to advect some of the particles. In the first particle tracking experiment, we advect the particles using the full, time-varying velocity fields, which have been saved every 5 model days. Using the full velocity fields to advect our Lagrangian particles provides us with a numerical estimate of the GLM velocity of the weakly stratified subtropical mode water. In the other, we use the Eulerian-mean velocity fields. Following Lagrangian particles in the Eulerian-mean velocity fields provides us with a velocity that is conceptually similar to the GLM velocity but with the effect of eddies removed. The difference between these two results is an estimate of the rectified effect of eddies on the bulk transport of fluid in the subtropical gyre.
The results of our Lagrangian particle tracking experiments are shown in Fig. 4. The particles advected by the full velocity fields remain above the thermocline, while the particles advected by the Eulerian-mean velocities are downwelled through the thermocline. The average vertical velocity of the particles in the Eulerian-mean velocity fields is initially comparable with the Ekman pumping velocity computed from the wind stress curl, but as they are downwelled the vertical motion slows to approximately half of the gyre-averaged Ekman pumping velocity.
We performed a similar experiment using a coarserresolution version of our model that included a Gent-McWilliams mesoscale eddy parameterization (Gent FIG. 4. The results of the particle tracking experiments, showing the location of the center of mass of a cloud of particles every 3 months. Purple circles show the center of mass of the particles in the full velocity field, and the green squares show the center of mass of the particles in the Eulerian-mean velocity field. The thin black line shows the trajectory given by the gyre average Ekman pumping velocity. et al. 1995). The model velocity fields showed downwelling, while the residual velocity fields, model plus bolus, did not. This result indicates that the GM mesoscale eddy parameterization is able to reproduce the eddy dynamics that are responsible for the cancellation of Ekman pumping.
The downwelling experienced by particles in the Eulerian-mean fields suggests that there is an Ekman cell in the Eulerian-mean velocities, while the trajectory of the particles in the full velocity fields supports the hypothesis that there exists an eddy-induced circulation that opposes the Eulerian-mean downwelling.
6. Vertical fluxes of potential vorticity Polton and Marshall (2007) use a framework that links vertical fluxes of potential vorticity with the advection of buoyancy to examine similarities between overturning circulations in the Southern Ocean and wind-driven gyres. Polton and Marshall (2007) find an eddy-mean balance in the Southern Ocean and in strongly eddying regions of the subtropical gyre but were hindered by the horizontal resolution of their model; at 1 /48, it was unable to satisfactorily resolve the mesoscale eddy field. Since our model is better able to resolve the mesoscale eddy field, we expect to find a greater contribution from eddy-induced fluxes of potential vorticity throughout the gyre.
Potential vorticity fluxes have been used extensively to investigate the structure of geophysical flows. Haynes and McIntyre (1990) laid the groundwork with their impermeability theorem by showing that potential vorticity substance is unable to pass through density surfaces. Schär (1993) generalized Bernoulli's theorem to show that in a statistically steady state, potential vorticity flux vectors are constrained, even in time-varying diabatic flows, to lie along the intersections between surfaces of constant density and surfaces of constant Bernoulli potential, here defined as P 5 (u 2 1 y 2 )/2 1 p/r 0 1 gz, where u and y are the horizontal velocity components, p is the pressure, r 0 is the reference density, and g is the gravitational constant. The vertical component of the velocity field is excluded to maintain consistency with the hydrostatic approximation. The impermeability of Bernoulli potential surfaces is less general than the impermeability theorem of Haynes and McIntyre (1990), since it applies only at a statistically steady state.
We use the potential vorticity flux framework described by Polton and Marshall (2007), based on earlier work by  and Polton and Marshall (2003), and the reader is directed to those papers or the attached appendix for a detailed discussion of the framework. This method essentially takes advantage of the impermeability of Bernoulli surfaces to potential vorticity flux vectors at a statistically steady state (Schär 1993) to define regions on a depth surface over which the area integral of vertical potential vorticity fluxes is constrained to sum to zero. While the net flux of potential vorticity is constrained to be zero, there is no other restriction on the relative magnitude of the different terms. We use this framework to explore the balance between Eulerian-mean advection and eddy-induced transport identified in section 5. In a statistically steady state, the vertical potential vorticity fluxes can be decomposed into the following different processes: Buoyancy changesJ buoy 5 q z B , Mechanical forcingJ fric 5 k 3 F Á =r, and (4) EddiesJ eddy 5 k 3 F R Á =r 2 q z = Á u 0 r 0 .
Here, the overbar represents Eulerian-mean terms, the tilde represents terms constructed from Eulerian-mean variables,Q 5 2(1/r 0 )q z (›r/›z) is potential vorticity, q z is the vertical component of Eulerian-mean absolute vorticity, w is the Eulerian-mean vertical velocity, B is the Eulerian mean of all nonadvective terms in the buoyancy equation, F is the Eulerian-mean momentum forcing, and F R is the Reynolds stress contribution to the Eulerian-mean momentum equation. These fluxes are integrated between closed contours of Eulerian-mean Bernoulli potential on a constant depth surface. Figure 5a shows a synopsis of the potential vorticity fluxes due to Eulerian-mean advection, eddies, and buoyancy forcing in the subtropical gyre. The negative fluxes from Eulerian-mean advection are consistent with Ekman-driven downwelling in the Eulerian-mean fields. The vertical potential vorticity flux intensity from eddies is shown in Fig. 5b and is predominantly balanced by fluxes due to Eulerian-mean advection, shown in Fig. 5c. The residual between eddies and Eulerian-mean advection is shown in Fig. 5d with a reduced color scale and is balanced by fluxes due to buoyancy forcing.
The potential vorticity fluxes in the subpolar gyre are shown in Fig. 6, with Figs. 6a and 6c providing a synopsis of the flux intensities plotted on two sets of axes for the sake of clarity. The positive values for Eulerian-mean advection are consistent with Ekman-driven upwelling. The Eulerian-mean advective fluxes, shown in Fig. 6d, are balanced by fluxes due to eddies (Fig. 6b) and buoyancy forcing.
Our results show that the potential vorticity flux due to time-mean advection is opposed by an eddy-induced flux and suggest that the rectified effect of eddies opposes the time-mean advection, consistent with the results of the Lagrangian particle tracking experiments.

Depth-density overturning streamfunctions
Our Lagrangian particle tracking experiments and analysis of potential vorticity fluxes have provided strong evidence for an eddy-driven circulation that opposes the Eulerian-mean Ekman downwelling. We now seek to quantify this effect through the use of an overturning streamfunction.
Initially, we sought to construct an overturning circulation around Eulerian-mean geostrophic streamlines, similar to the method described by Marshall and Radko (2003) for the Antarctic Circumpolar Current. However, both the shape and topology of the gyres in our model, and in the ocean, vary substantially with depth, complicating the definition of an unambiguous coordinate system to serve as the gyre equivalent of the depth-latitude coordinate system used for the Antarctic Circumpolar Current. Instead, we use the depth-density overturning streamfunction developed by Nurser and Lee (2004). Using the depth-density formalism obviates the need to specify an explicit geometry for our streamfunction. Since our model employs a linear equation of state that depends only on temperature, there is a one-to-one correspondence between density and temperature. Henceforth, we will refer to these as depth-temperature streamfunctions.
The depth-temperature overturning streamfunction calculated using the Eulerian-mean velocity and temperature fields is shown in Fig. 7a. There is evidence of two circulations superimposed on top of each other. The blue lobes, rotating counterclockwise, represent a meridional overturning circulation driven by dense-water formation in the north of the domain and upwelling in the sponge at the southern boundary along the imposed temperature profile shown by the blue line. This meridional overturning circulation spans the entire domain and will henceforth be referred to as the MOC. The red lobes, which rotate clockwise, show fluid sinking in the 158-208C range and upwelling at colder temperatures. When we calculate this streamfunction using the instantaneous fields, and then average, we obtain a strikingly different result. The 12-Sv (1 Sv [ 10 6 m 3 s 21 ) MOC in Fig. 7b is stronger than the circulation of 6 Sv from Fig. 7a. The overturning streamfunction calculated from the instantaneous fields is also more adiabatic; the water mass transformations occur primarily at the surface or in the southern sponge region.
Since the only difference between Figs. 7a and 7b is the time-varying component, this provides a way to quantify the rectified transport due to eddies. We define our eddydriven depth-temperature overturning streamfunction as where the overbar again represents a time-averaged operator, and C inst is calculated at each time level in the saved data and then averaged. The eddy-driven overturning streamfunction is shown in Fig. 7c. The Ekman cell in Fig. 7a is somewhat obscured by the MOC. We can remove the MOC by appropriately masking the domain when calculating the depthtemperature streamfunction. Masking the sponge removes the upwelling component of the MOC. It is less straightforward to mask the downwelling limb of the MOC. We have tested both latitude-based criteria, masking northward of a specified latitude, and temperature-based criteria, masking regions with a temperature that is less than a specified cutoff temperature. Neither of these methods was able to completely isolate the Ekman cell from the MOC, but both gave similar results. Figure 7d shows the depth-temperature overturning streamfunction FIG. 7. Depth-temperature overturning streamfunctions from our idealized model. The black arrows indicate the direction of fluid transport: red lobes rotate clockwise, and blue ones rotate counterclockwise. The blue line represents the temperature profile in the sponge region and is thus the profile along which we expect upwelling in the sponge region. (a) Calculated from the Eulerian-mean fields and shows the meridional overturning circulation in blue with a second circulation, the Ekman cell, superimposed in red. (b) The depth-temperature streamfunction calculated from instantaneous model fields and then averaged. (c) The difference between these two is the eddy contribution to the depthtemperature streamfunction. calculated excluding the southern-and northernmost 108 of latitude. In the masked streamfunction, the two red lobes from Fig. 7a are now connected and form a single, thermally indirect circulation. The masked streamfunction does not close because we are considering an open domain. We have truncated the streamfunction at the dashed line in Fig. 7d, which represents the maximum temperature in our model at each depth. Figure 7e shows the masked instantaneous depthtemperature overturning streamfunction. The sinking at approximately 158C occurs within the subtropical gyre and is consistent with the expected 10 Sv of Ekman pumping from our wind stress curl. However, as shown by this streamfunction, the Ekman pumping is unable to penetrate deeper than 200 m below the surface. The other feature to note is that this circulation is thermally direct; cold water is downwelled and warmer water is upwelled. In this framework, the eddy-induced overturning is stronger than the Eulerian-mean overturning, leading to the subduction of colder water and the upwelling of warmer water. Figure 7f shows the eddyinduced overturning streamfunction calculated in the masked domain.
The plots shown in Fig. 7 exclude the deep ocean below 600-m depth. The circulation below this depth is less interesting, with the unmasked streamfunctions (Figs. 7a-c) closing between 1000-and 2000-m depth. The masked streamfunctions (Figs. 7d and 7e) continue to similar depths but do not close because of the masking procedure, while Fig. 7f closes in the upper 600 m.

Concluding remarks
We found three lines of evidence to support our hypothesis that an eddy-driven circulation opposes the Eulerian-mean Ekman cell: Our Lagrangian particle tracking experiments show that the residual vertical velocity averaged over the subtropical gyre is small; the potential vorticity flux budget shows a balance between Eulerian-mean advection and eddies, with a small, but nonnegligible, contribution from buoyancy changes; and in depth-temperature space, overturning streamfunctions indicate that an eddy-driven circulation opposes the Eulerian-mean Ekman cell and prevents the downwelled fluid from penetrating deeper than 200 m below the surface.
The main conclusion of this paper is that eddies play a substantial role in setting the residual overturning circulation within wind-driven gyres, which have long been described with time-mean linear theories. While the idealized nature of our model may have enhanced this effect and contributed to the near-complete cancellation of the Ekman-driven downwelling, we expect some level of cancellation to be a robust result. Near-complete cancellation between eddies and the Eulerian-mean flow is consistent with the results of Griffies et al. (2015), who found a balance between vertical heat fluxes from Eulerian-mean advection and eddy advection throughout their eddy-resolving global model.
The presence of a nonnegligible contribution from buoyancy changes in our potential vorticity flux budget is consistent with the results of Henning and Vallis (2004), who found a partial balance between diapycnal eddy-driven mass fluxes and mean upwelling. They also found that the residual from this balance was smaller than either of the terms and was balanced by diffusion. Our results are also consistent with the idealized numerical experiments of Radko and Marshall (2004a), who found a balance between mean advection and eddies in global buoyancy and potential vorticity budgets.
Although we show compensation between an eddydriven overturning circulation and an Eulerian-mean circulation, we wish to highlight that this effect is subtly different to the concept of ''eddy compensation.'' Within the Southern Ocean literature, the term eddy compensation is used to refer to the dynamic change of the eddy-driven overturning circulation in response to a change in the wind stress (Viebahn and Eden 2010). The eddy-driven overturning circulation partially compensates for changes in the surface Ekman transport (Morrison and Hogg 2013). Eddy compensation is part of the dynamic adjustment of the residual overturning circulation toward a new equilibrium. Unlike eddy compensation, the phenomenon discussed in this paper is an equilibrium balance between an Eulerian-mean circulation and an eddy-induced circulation. To emphasize the difference between the dynamic concept of eddy compensation and the equilibrium concept discussed here, we suggest eddy cancellation as the name for the phenomenon discussed in this paper.
Sverdrup balance predicts that the horizontal gyre transport will increase with the wind stress curl, but it is unclear how the vertical transport would respond or the time scale over which this response would occur. Within the Southern Ocean literature, there are numerous studies showing that the presence of eddies in models of the Southern Ocean reduces the sensitivity of the residual meridional overturning circulation to changes in wind stress (see, e.g., Abernathey et al. 2011;Hallberg and Gnanadesikan 2006;Henning and Vallis 2005;Munday et al. 2013;Viebahn and Eden 2010). It may therefore be of interest to investigate how the Eulerian-mean downwelling and eddy-driven upwelling in eddy-resolving models of subtropical gyres respond to changes in the wind stress. In particular, understanding the time scale of the eddy-driven response should provide insight into the dynamical and biogeochemical response of subtropical gyres to changes in the wind stress.
The existence of an eddy-driven circulation that opposes the Eulerian-mean Ekman cell would explain the unexpected eddy-driven downwelling in the subpolar gyre found by McGillicuddy et al. (2003). Since our results show that an eddy-driven circulation opposes the Ekman velocity, we expect to see eddy-driven downwelling in the subpolar gyre. Our result may also help to explain the low correlation between maps of biological activity and time-mean Ekman pumping (see, e.g., Palter et al. 2005).
If large-scale, Ekman-driven downwelling is effectively canceled at depth in subtropical gyres, this has substantial implications for our understanding of the nutrient budget of these regions. Ekman-driven downwelling is invoked as a mechanism to remove nutrients from the subtropical gyres (e.g., Jenkins and Doney 2003;Kähler et al. 2009;McGillicuddy et al. 1998;Williams and Follows 1998). The presence of a gyre-scale rectified transport opposing this downwelling would allow nutrients to be recycled more easily within subtropical gyres, thereby reducing the quantity of new nutrients required to supply the observed levels of productivity. Although Lee and Williams (2000) examined the balance between eddy diffusion, eddy advection, and Ekman transport in the horizontal movement of nutrients, the idea of eddy cancellation in the vertical is a qualitatively different role for eddies than has been previously proposed (e.g., McGillicuddy et al. 1998McGillicuddy et al. , 2003McGillicuddy 2016). By lowering the required supply of nutrients, this could settle the longstanding discrepancy between estimates of nutrient supply and export production in subtropical gyres. Results from such a study will be reported in a future manuscript. Future Fellowship FT120100842.
Configuration and input files for our model as well as copies of the analysis scripts can be found online (at http://doddridge.me/publications/dmh2016/).

Calculating the Potential Vorticity Fluxes
Following  and Polton and Marshall (2007), we calculate vertical fluxes of potential vorticity within closed contours of Eulerian-mean Bernoulli potential.  constructed an integral constraint on vertical fluxes of potential vorticity within closed contours of Bernoulli potential. The framework developed by  was used by Polton and Marshall (2007) to explore the vertical structure of subtropical gyres and to examine dynamical similarities between these gyres and the Southern Ocean.
Under the Boussinesq approximation with an incompressible fluid in which the density r is a materially conserved tracer, the potential vorticity evolution equation from Ertel (1942) becomes where r 0 is a reference density; Q 5 2(q/r 0 ) Á =r is the potential vorticity; q is the absolute vorticity; B represents all nonadvective processes that affect density; F represents the mechanical and body forcings acting on the fluid, such as viscosity and surface wind forcing; and g is the gravitational acceleration. Because our simulations employ a single variable linear equation of state, we use in situ density for these calculations. In models with more complex thermodynamics, or in the ocean, potential density or neutral density (McDougall 1987) would be more appropriate.
The negative sign in the definition of potential vorticity ensures that the large-scale potential vorticity is positive in the Northern Hemisphere and hence that the advective flux in the simulations we consider is in the same direction as the fluid velocity, rather than antiparallel. Throughout this paper the potential vorticity Q is approximated as where f and z represent the vertical components of planetary and relative vorticity, respectively. For particularly large-scale or weak-flow applications, this may be simplified even further by ignoring the contribution from relative vorticity. As noted by Haynes andMcIntyre (1987, 1990), Eq. (A1) may be recast into a flux form equation that relates the local time derivative of potential vorticity to the divergence of a potential vorticity flux vector J.
The flux form of the potential vorticity conservation equation is where J is the potential vorticity flux vector that includes advective and nonadvective flux terms. As defined in Eq. (A3), only the divergence of the potential vorticity flux vector is constrained. Therefore, the potential vorticity flux may also contain any arbitrary divergence free vector field; this represents a gauge invariance for the potential vorticity flux vector field. Since the divergence of the curl of an arbitrary vector field is identically zero, we include a gauge term in this form. The resulting potential vorticity flux vector field is given by where = 3 A represents the gauge. The form presented in Eq. (A4) explicitly shows the advective and nonadvective contributions to the potential vorticity flux.
To derive an alternate formalism, we begin with the hydrostatic Boussinesq momentum equation written in terms of the Bernoulli potential and take the cross product of each term with the gradient of density. After some simplification, this gives where P 5 (u 2 1 y 2 )/2 1 p/r 0 1 gz is the Bernoulli potential, consistent with the hydrostatic approximation and v 5 (u, y, 0). This form makes explicit the impermeability of density and Bernoulli surfaces in the steady-state solution, provided that the gauge is chosen such that = 3 A is equal to zero. Haynes andMcIntyre (1987, 1990) were the first to generalize Ertel's theorem to show that isentropic surfaces are impermeable to potential vorticity. This is known as ''the impermeability theorem.' ' Schär (1993) generalized Bernoulli's theorem to show that in a statistically steady state, potential vorticity flux vectors are constrained to lie along the intersections between surfaces of constant density and surfaces of constant Bernoulli potential, even in timevarying diabatic flows. The impermeability of Bernoulli potential surfaces is less general than the impermeability theorem of Haynes and McIntyre (1990), since it applies only at a statistically steady state.

a. Eulerian-mean potential vorticity flux
Equations (A4) and (A5) describe the instantaneous flux of instantaneous potential vorticity. It is possible to examine Eulerian-mean potential vorticity fluxes by applying a suitable time filter to the instantaneous equations. To derive the Eulerian-mean potential vorticity fluxes, we start by filtering the momentum and buoyancy equations into time-mean X and timevarying X 0 components. Using this notation, the timefiltered hydrostatic Boussinesq momentum equation takes the form whereP represents the Bernoulli potential constructed from Eulerian-mean quantities, andP 5 (u 2 1 y 2 )/2 1 p/r 0 1 gz, where u and y are the Eulerian-mean horizontal velocity components, p is the Eulerian-mean pressure, r 0 is the reference density, g is the gravitational constant, z is the vertical coordinate, and the time derivative represents changes over time scales comparable with, or longer than, the time scale of the time-mean operator. Equation (A6) can be simplified if the eddy contributions are grouped into a single term Applying the same time filter to the buoyancy equation results in where the nondivergence of the velocity field has been used to simplify the second term on the righthand side. If these time-filtered momentum and buoyancy equations are used to construct the potential vorticity flux vectors, then we get These equations have the same form as the instantaneous equations, shown in Eqs. (A4) and (A5), except that F and B have been replaced by F 1 F R and B 2 = Á u 0 r 0 , respectively, to take into account the rectified effect of eddies. Since they are exact, these potential vorticity fluxes are necessarily equivalent to the results derived in Young (2012) and Maddison and Marshall (2013).
To preserve the impermeability of surfaces of constant density and Bernoulli potential, the gauge used throughout this paper is chosen such that = 3 A 5 0.

b. Vertical PV flux
Following , we use the impermeability of surfaces of constant density and Bernoulli potential to create volumes into which there is no net flux of potential vorticity. We construct these volumes by taking an isosurface of Bernoulli potential and enclosing a volume with a surface of constant depth. Since there can be no net flux of potential vorticity across the Bernoulli potential isosurface, there must also be no net flux across the depth surface that forms the closed volume, provided that the amount of potential vorticity substance inside the volume is constant. While the net flux is constrained to zero by construction, it is not clear a priori which processes should dominate the dynamical balance.
The integral may be performed with either of the expressions for potential vorticity flux, Eqs. (A9) and (A10), or any linear combination of the two.
When performing this integral, the =P 3 =r term vanishes. This was shown by , who, through the use of vector identities and Stokes' theorem, converted an area integral inside a closed contour of Bernoulli potential, ðð l (=P 3 =r) dA, into a line integral around contours of constant Bernoulli potential with the form 2 þ ›l r=P Á dr , which must equal zero, since the edge of the area is defined as a line of constant Bernoulli potential, and so =P Á dr is identically zero everywhere on ›l. The only other terms in Eq. (A10) include time derivatives, which can be used to check the impact of model drift on the results. These drift terms relate to changes in the storage of potential vorticity substance within the constructed volume. We obtain our potential vorticity flux diagnostics by integrating the vertical component of Eq. (A9) minus Eq. (A10) between closed contours of Bernoulli potential on a depth surface. The integral has the form ðð lk Á r 0Q u 1 q(B 2 = Á u 0 r 0 ) 1 F 1 F R 2 Dr r 0 gk 3 =r 2 q ›r ›t 2 ›v ›t 3 =r dA 5 0, where l is defined as the area between two closed contours of Bernoulli potential. The terms in Eq. (A13) describe the vertical potential vorticity flux due to specific physical processes. This can be made explicit by rewriting Eq. (A13) as ðð l (J adv 1J buoy 1J fric 1J eddy 1J drift ) dA 5 0, whereJ adv 5 r 0Q w (A15) is the potential vorticity flux from advection by the Eulerian-mean vertical velocity, is the flux due to buoyancy forcing, is the flux due to forcings in the momentum equation, J eddy 5 k 3 F R Á =r 2 q z = Á u 0 r 0 (A18) is the eddy-induced flux, and J drift 5 2k 3 ›v ›t Á =r 2 q z ›r ›t (A19) is the flux due to model drift. In each of these, the superscript z refers to the vertical component.

c. Evaluating potential vorticity fluxes in the model
The fluxes shown in Eqs. (A15)-(A19) were evaluated from a simulation spanning 10 yr of model time. The potential vorticity and Bernoulli potential were constructed from the Eulerian-mean variables. Consistent with the hydrostatic approximation, the vertical velocity was not included in the calculation of the Bernoulli potential. The eddy contribution to the potential vorticity flux J eddy 5 k 3 F R Á =r 2 q z = Á u 0 r 0 (A20) consists of the rectified effect of eddies in the momentum equation and the buoyancy equation. These rectified eddy terms were calculated online by averaging and saving the terms q z u, q z y, uu, yy, and ur. The eddy contributions were then calculated as the difference between these terms and the product of the individual Eulerian-mean fields. For example, the zonal velocity eddy term is calculated as u 0 u 0 5 uu 2 u u , which requires that we average uu and u and save them at the end of the 10-yr simulation. Unlike Polton and Marshall (2007), we use the full three-dimensional velocity to calculate the divergence of the eddy buoyancy flux, since we were able to save these variables. The potential vorticity fluxes due to model driftJ drift were evaluated and consistently found to be negligible at the depths analyzed in this paper.