## Abstract

An eddy-resolving numerical model of a zonal flow, meant to resemble the Antarctic Circumpolar Current, is described and analyzed using the framework of J. Marshall and T. Radko. In addition to wind and buoyancy forcing at the surface, the model contains a sponge layer at the northern boundary that permits a residual meridional overturning circulation (MOC) to exist at depth. The strength of the residual MOC is diagnosed for different strengths of surface wind stress. It is found that the eddy circulation largely compensates for the changes in Ekman circulation. The extent of the compensation and thus the sensitivity of the MOC to the winds depend on the surface boundary condition. A fixed-heat-flux surface boundary severely limits the ability of the MOC to change. An interactive heat flux leads to greater sensitivity. To explain the MOC sensitivity to the wind strength under the interactive heat flux, transformed Eulerian-mean theory is applied, in which the eddy diffusivity plays a central role in determining the eddy response. A scaling theory for the eddy diffusivity, based on the mechanical energy balance, is developed and tested; the average magnitude of the diffusivity is found to be proportional to the square root of the wind stress. The MOC sensitivity to the winds based on this scaling is compared with the true sensitivity diagnosed from the experiments.

## 1. Introduction

Changes in wind stress over the Southern Ocean may be responsible for modulating the strength of the global meridional overturning circulation (MOC) (Toggweiler 2009). Such wind-induced changes in the MOC could help regulate glacial–interglacial cycles by venting CO_{2} from the deep ocean to the atmosphere (Toggweiler and Russell 2008; Anderson et al. 2009; Marshall and Speer 2011). The mechanism could also play an important role in future climate change; the westerlies appear to be shifting south because of greenhouse gas emissions and ozone depletion (Thompson and Wallace 2000; Marshall 2003; Polvani et al. 2011), and Toggweiler and Russell (2008) hypothesize that in response the MOC will strengthen, but by how much? To what extent is the Southern Ocean MOC controlled by the winds?

Since Johnson and Bryden (1989), we have recognized the existence of an eddy-driven overturning circulation in the Southern Ocean potentially large enough to completely cancel the wind-driven Ekman overturning. The actual MOC is the small residual between these two opposing circulations. Work by Toggweiler and Samuels (1998), Speer et al. (2000), and Marshall and Radko (2003, hereafter MR03) showed that, for realistically weak values of interior diapycnal mixing, the residual overturning transport in the subsurface Southern Ocean must proceed along mean isopycnal surfaces. The residual circulation can cross isopycnals in the surface diabatic layer, where cross-isopycnal advection can be balanced by direct diabatic forcing from the atmosphere (Marshall 1997). Therefore, from a diagnostic point of view, the strength and sense of the MOC can be inferred from surface buoyancy-flux data, as done by Speer et al. (2000) and Karsten and Marshall (2002b), independently of the wind stress. This thermodynamic perspective also implies that the MOC is sensitive to surface buoyancy fluxes, as hypothesized by Watson and Naveira Garabato (2006) or Badin and Williams (2010). Our goal here is to study the relationship between wind stress, overturning circulation, and surface buoyancy flux in a model that explicitly resolves mesoscale eddies, bypassing the need for any a priori assumptions about the eddy response.

On a related note, it is well established that coarse-resolution ocean models do not accurately simulate the response of the Southern Ocean overturning to changes in wind stress forcing when compared with eddy-resolving models. This is true of both realistic models (Hallberg and Gnanadesikan 2006; Farneti et al. 2010) and models with simplified geometry and forcing (Henning and Vallis 2005). In general, models that permit eddies seem less sensitive to changes in wind, whether the focus is the overturning circulation (as in the above works), the zonal transport (Hutchinson et al. 2010), or the transport of tracers such as anthropogenic carbon (D. Munday 2011, unpublished manuscript). Most of these results are ultimately due to compensation between the wind- and eddy-driven overturning circulations, which is more complete when mesoscale eddies are explicitly resolved rather than parameterized. The lack of a robust parameterization for mesoscale eddies is indicative of our incomplete understanding of the nature of eddy-driven circulations. Most recently, Viebahn and Eden (2010) studied the sensitivity of the residual MOC to the wind in an idealized model and found that changes in eddy kinetic energy (EKE) and eddy diffusivity play a central role in determining how the compensation occurs.

Our goal in this study is to further explore the physical mechanisms that determine the sensitivity of the residual MOC to changes in wind forcing. In particular, there are two questions not previously addressed that we wish to pursue here. First is the influence of the boundary condition for buoyancy. Second, we wish to develop a simple theory based on physical principles capable of explaining the MOC sensitivity. To study these issues, we reduce the system to its essential elements: an Ekman-driven and an eddy-driven circulation in a zonally symmetric channel with buoyancy forcing at the surface. This system was studied analytically by MR03, who invoked a closure for the eddies, but here we realize it as a high-resolution numerical model. The strength of the Ekman circulation obviously depends linearly on the winds; the strength of the eddy-driven circulation is determined by the geostrophic turbulent dynamics of the model. We vary the strength of the wind stress and diagnose the steady-state residual overturning circulation.

We find that increased eddy circulation does generally compensate for increased Ekman circulation under stronger winds. However, the degree of compensation depends on the surface boundary conditions. When the surface heat fluxes are held fixed, the residual MOC strength is relatively insensitive to the winds. With an interactive heat flux, we recover the results of Viebahn and Eden (2010): a residual MOC that increases weakly with the winds and whose sensitivity is set primarily by changes in eddy diffusivity. We develop a scaling theory for the eddy diffusivity dependence on the wind and apply this scaling to reconstruct the eddy response. This method yields a closed theory for the sensitivity of the residual MOC, which, despite many approximations, shows encouraging agreement with the results from the numerical model.

Section 2 describes the model setup, a reference solution, and the basic experimental results under differing values of wind stress. In section 3, we analyze the results in terms of the buoyancy budget and discuss the constraints imposed by the surface boundary condition for buoyancy. Section 4 describes a framework for understanding the MOC changes in terms of changes in Ekman circulation, isopycnal slope, and eddy diffusivity. Our scaling for the eddy diffusivity and the resulting MOC sensitivity estimates are then presented. We summarize the results and discuss their connection with the real ocean in section 5.

## 2. Experiments with numerical model

### a. Modeling philosophy

The Southern Ocean is dominated by the Antarctic Circumpolar Current (ACC), a strong eastward flow in thermal-wind balance with the strong density front separating polar from tropical waters (Rintoul et al. 2001). This flow circumnavigates the globe and connects back on itself, inspiring a comparison with the large-scale atmospheric jets (Thompson 2008). Strong atmospheric westerly winds blow over the surface, driving an equatorward Ekman flow. The surface buoyancy flux—a combination of radiative, latent, and sensible heat fluxes as well as freshwater fluxes from evaporation, precipitation, and ice-related processes—is notoriously uncertain because of poor data sampling (Cerovečki et al. 2011). Nevertheless, the general pattern (shown in Fig. 1) indicates buoyancy loss in the extreme south polar regions, buoyancy gain on the poleward flank of the ACC, and buoyancy loss in some regions on the equatorward flank associated with mode water formation. Although the current meanders and splits as it makes its way around topographic features, authors such as de Szoeke and Levine (1981) and Ivchenko et al. (1996) have argued that, when the real ACC is described using a “streamwise average” view, the large-scale dynamics bear a close resemblance to zonally symmetric models.

Indeed, zonal channel models with highly idealized geometry form the foundation of contemporary theories of the Southern Ocean circulation, capturing the essential physics of the system and providing insight into important mechanisms (Munk and Palmén 1951; McWilliams et al. 1978; Marshall 1981; Johnson and Bryden 1989; Marshall 1997; Olbers et al. 2004; Marshall and Radko 2006, among many). The Southern Ocean MOC, however, exports and imports water from other ocean basins (e.g., Ganachaud and Wunsch 2000; Talley 2008). Antarctic Bottom Water (AABW) flows out of the Southern Ocean in the deepest layers. North Atlantic Deep Water (NADW) and Circumpolar Deep Water (CDW) flow in (poleward) at intermediate depths, and Antarctic Intermediate Water (AAIW) and Subantarctic Mode Water (SAMW) flow equatorward in the upper thermocline. As a result, channel-only models that attempt to investigate the Southern Ocean MOC without representing other basins find vanishingly weak deep residual circulations (Karsten et al. 2002; Kuo et al. 2005; Cessi et al. 2006; Cerovečki et al. 2009). Some authors have tackled this problem by attaching closed basins to their channels. This approach can certainly yield insights, but it also adds to the complexity of the problem by introducing gyre dynamics. When such basins are global scale, as in Wolfe and Cessi (2009), the computational cost of an eddy-resolving model becomes immense. When they are small (on the same order of the channel itself), as in Henning and Vallis (2005) and Viebahn and Eden (2010), the link with the real ocean is less clear.

We choose to address this problem in a novel way: by including a narrow “sponge layer” along the channel’s northern boundary, in which the temperature is relaxed to a prescribed exponential stratification profile. This diabatic forcing provides a return pathway for deep residual overturning, which otherwise would not be able to cross isopycnals. Physically, the sponge layer encapsulates all the diabatic processes occurring outside of the Southern Ocean, such as deep-water formation by air–sea heat fluxes in the North Atlantic or diapycnal mixing in the abyss. The disadvantage of this method is that the stratification at the northern boundary cannot change significantly. The advantage is that it provides a clean, simple framework in which to investigate nonzero residual circulations, focusing on the dynamics of the channel alone rather than the complex teleconnections of the global problem (Wolfe and Cessi 2011). In combination with appropriate surface wind and buoyancy forcing, we will see that this configuration can produce realistic overturning cells.

Given the many idealizations made in constructing our model, we must interpret our results with care. We emphasize that our goal is not to make quantitative predications for the real global ocean–atmosphere system; rather, we hope to gain insight into the underlying physical mechanisms that govern this system in order to inform the interpretation of more realistic models and observations.

### b. Model physics and numerics

The basic physical system simulated by our model is a Boussinesq fluid on a beta plane with a linear equation of state and no salinity. The model is forced mechanically by a surface stress and thermodynamically by a surface heat flux as well as by the aforementioned sponge-layer restoring. Mechanical damping is provided by linear bottom drag; there is no topography. Key physical and numerical parameters are given in Table 1.

The surface thermal forcing in our model is intended to mimic, in a simplified way, the observed buoyancy flux over the Southern Ocean (see Fig. 1). In the first set of experiments, a heat flux is simply prescribed to include a region of cooling in the far south of the domain, heating in the middle, and cooling again farther north. These regions are intended to represent the buoyancy loss associated with AABW formation, the buoyancy gain over the ACC, and the buoyancy loss associated with AAIW–SAMW formation north of the front, respectively. More precisely, the heat flux has the form

and *Q* = 0 north of this point, with *Q*_{0} = 10 W m^{−2}. The term *L _{y}* is the length of the channel in

*y*(

*Q*is positive downward: i.e., heat flux into the ocean). This simple pattern of buoyancy flux is consistent with a recent review of all the available air–sea buoyancy flux data products by I. Cerovečki et al. (2011, personal communication).

Inside the sponge layer, the temperature *T* is relaxed to the prescribed temperature profile

which describes an exponential decay from Δ*T* at the surface to 0 at depth −*H* with a scale height of *h*. The relaxation coefficient increases from 0 (meaning no relaxation) at the southern edge of the sponge layer (*y* = *L*_{sponge}) to 7 day^{−1} at the northern boundary (*y* = *L _{y}*). The choice of an exponential temperature profile was motivated by observations (Karsten and Marshall 2002a), laboratory studies (Cenedese et al. 2004), and modeling results (Karsten et al. 2002; Henning and Vallis 2005; Wolfe and Cessi 2009). The results described in the rest of the paper all use

*h*= 1000 m, a value close to the “natural” stratification that arises when the sponge layer is turned off and to the observed stratification on the equatorward flank of the real ACC. We experimented with several values of the stratification depth

*h*and found that the MOC transport was rather insensitive to this choice.

The final key element of the forcing is the wind stress. A zonal stress is applied at the surface of the form

For the base-case simulation, *τ*_{0} = 0.2 N m^{−2}, but a central point of our study is to explore the strength of the MOC given different values of *τ*_{0}.

Dissipation is mainly accomplished through linear bottom drag. A stress is applied at the bottom of the form

where *r _{b}* is a bottom drag coefficient and

**u**

*the horizontal component of the bottom velocity.*

_{b}The model code is the Massachusetts Institute of Technology general circulation model (MITgcm), a general-purpose primitive equation solver (Marshall et al. 1997a,b). The domain is a Cartesian grid 1000 km long (i.e., zonal direction), 2000 km wide (i.e., meridional direction), and 2985 m deep. Although this domain is relatively narrow, the zonal symmetry means that a larger domain would not alter the results and would only add computational cost. The domain size does not appear to constrain the eddy size, because a typical eddy size is ~200 km. We resolve the first baroclinic deformation radius (approximately 15 km in the center of the domain), employing 5-km horizontal resolution and with 30 vertical levels, with spacing increasing from 10 m at the surface to 280 m at the bottom. A realistically effective diapycnal diffusivity (*κ _{υ}* = 0.5 × 10

^{−5}m s

^{−2}) is maintained thanks to the second-order-moment advection scheme of Prather (1986) (see also Hill et al. 2011). To maintain a surface mixed layer, we employed the

*K*-profile parameterization (KPP) mixing scheme (Large et al. 1994). In our case, this scheme simply acts to mix tracers and momentum over a layer of roughly 50-m depth.

The model was spun up from rest for approximately 200 yr until it reached a statistically steady state, as indicated by the mean kinetic energy. A typical eddy temperature field from the equilibrated state is shown in Fig. 2. Averages were performed over 20-yr intervals. In cases where parameters were changed, the model was allowed to reach a new equilibrium before taking an average.

### c. The zonal momentum balance

Because the meridional flux of momentum by Reynolds stresses is relatively small, the depth-integrated zonal-average momentum balance dictates (cf. Cessi et al. 2006) that

where the overbar indicates a zonal and time average. This balance states that the momentum input by the wind (constant in *x* and time) is balanced by bottom drag on a mean zonal flow at the bottom. In the real ocean, in contrast, topographic form drag is believed to balance the wind stress (Munk and Palmén 1951; Johnson and Bryden 1989; Hughes 1997; Olbers 1998; Ferreira et al. 2005). This means that our model requires a significant steady bottom flow (~17 cm s^{−1}) and thus has an unrealistically large zonal transport: 788 Sv (1 Sv ≡ 10^{6} m^{3} s^{−1}) for the reference case. But most of this transport is barotropic and simply translates the entire system westward without any consequences for the overturning circulation. The zonal transport by the baroclinic flow is only 99 Sv.

A steady meridional circulation exists in Coriolis balance with these steady zonal stresses. Outside of the Ekman layers, this circulation is described by the streamfunction

where and . The absence of topography means that the surface Ekman flow is returned in a bottom Ekman layer, rather than by a geostrophic flow below topography. However, the strength of is independent of the nature of the bottom drag and is driven solely by the wind.

Likewise, as discussed in detail in section 4, the barotropic component of the flow does not participate in the eddy energy cycle, and thus we expect the eddy-driven circulation to be similar with or without topography. Experiments performed with a topographic ridge (but not described further here) support the conclusion that the presence of topography strongly damps the barotropic zonal flow but does not affect the MOC. We therefore expect that conclusions drawn from our model about the MOC can still apply to the real Southern Ocean, especially to the portion of the flow that occurs above major topographic features.

### d. Residual overturning circulation

To diagnose the residual MOC, we computed a streamfunction from the time- and zonal-mean transport in isopycnal layers, defined as

where *h* = −∂*z*/∂*b* is the layer thickness and *b*′ is a dummy variable of integration. (In practice, the average was performed in 22 discrete, uniformly spaced temperature layers.) This technique has become widely used for diagnosing transport in the presence of eddies (Döös and Webb 1994; Henning and Vallis 2005; Hallberg and Gnanadesikan 2006; Wolfe and Cessi 2009, 2010). The transport thus computed includes both the Eulerian-mean (Ekman) transport and the eddy-driven component. We can map this streamfunction into *z* coordinates by using the mean depth of buoyancy surfaces, . The leading-order equivalence between Ψ_{iso} in *z* coordinates and the transformed Eulerian-mean (TEM) residual circulation is well documented (Andrews et al. 1987; McIntosh and McDougall 1996; see section 4 for more on TEM theory). The most climate-relevant quantity is Ψ_{iso}, because it describes the circulation that advects tracers such as heat and carbon. Henceforth, when we refer to the MOC, we will generally be talking about Ψ_{iso} as defined in (7).

The MOC is characterized by three distinct cells, as shown in Fig. 3. In the interior of the domain, away from the surface and the sponge layer, the MOC is directed along mean isopycnals: that is, . Although the circulation is highly idealized, it shares several important features with the real Southern Ocean MOC, as described, for example, by Rintoul et al. (2001), Lumpkin and Speer (2007), Talley (2008), or Marshall and Speer (2011). The magnitude of Ψ_{iso} (~0.5 Sv) is realistic: if our channel were as long as the real Southern Ocean (a factor of about 25), the transport would be roughly 12 Sv. The broad upwelling band at middepth can be thought of as NADW–CDW. This upwelling water splits into two separate cells. The upper branch travels north, eventually encountering a region of cooling and subduction. The subduction in this northern region along the 4°C isotherm, driven by surface heat loss and accompanied by low values of Ertel potential vorticity, is reminiscent of SAMW–AAIW formation (McCartney 1977; McCarthy and Talley 1999). The water associated with downwelling in the far south of the domain resembles AABW in some respects; it is formed by buoyancy loss and is the coldest, densest water in the model. Given the complex physics of AABW formation on the continental shelf, the fact that much of the AABW circulates at depths blocked by topography, and the importance of diapycnal mixing for the lower limb overturning (Ito and Marshall 2008), this lower cell is not meant to be a truly realistic representation of AABW. All the overturning cells have an adiabatic pathway in the ocean interior and close diabatically in the sponge layer.

The surface heat flux is specified as a fixed function of latitude; consequently, the heat flux is felt by all isopycnals that graze the surface at that particular latitude. The cumulative distribution function (CDF) of surface temperature tells how likely a particular temperature is to be found at the surface and thus be exposed to diabatic transformation. Superimposed on Fig. 3 are the 5% and 95% values of *T* from the surface temperature CDF. (The mean SST is very close to the median value.) Nearly all of the diabatic MOC transport (i.e., advection across mean isopycnals) takes places in between these values. When plotted in *z* coordinates, the 95% CDF value is an effective measurement of the depth of the surface diabatic layer; below it, the contours of Ψ_{iso} and coincide. This view also reveals that the northernmost, shallow, counterclockwise MOC cell is contained almost entirely by the diabatic layer. We do not focus further on this shallow cell, concentrating from now on only on the two cells (lower and upper) that enter the adiabatic interior.

This definition of the MOC streamfunction [Eq. (7)] should be distinguished from the steady, Eulerian-mean overturning streamfunction [Eq. (6)]. The difference between the two circulations we define as the eddy circulation,

### e. Model response to wind changes

We now examine the MOC sensitivity to altered wind stress. We consider two principal cases. First, the surface buoyancy flux is held fixed as the winds are varied. In the second set of experiments, we employ an interactive, relaxation-type boundary condition. The fixed-flux boundary condition is a justifiable one for freshwater and incoming shortwave radiation, but an interactive boundary condition is more appropriate for sensible and latent heat (Haney 1971). The results are summarized in Fig. 5, where the strengths of the upper and lower cells in each experiment are plotted on a single graph. Because Ψ_{iso} is roughly constant along isopycnals below the diabatic layer, we diagnosed the transports by simply finding the maximum and minimum values of Ψ_{iso} below 500 m at *y* = 1800 km, 100 km south of the edge of the sponge layer. We will henceforth refer to these maximum and minimum values of Ψ_{iso} as MOC_{upper} and MOC_{lower}. Besides the individual upper and lower cells, there is a third relevant quantity: the total volume flux of upwelled deep water, MOC_{upwell} = MOC_{upper} − MOC_{lower}. This value is also shown in Fig. 5, along with the strength of , the Ekman circulation. In almost all cases, MOC_{upper}, MOC_{lower}, and MOC_{upwell} are weaker than .

In general, the various MOC values appear to have linear dependence on the wind. This is not a universal rule for all possible models and ranges of parameters (e.g., Viebahn and Eden 2010), but it is an accurate and useful approximation for our particular experiments. This simplification allows us to characterize the MOC sensitivities in a single number by a simple least squares linear fit applied to Fig. 5. The slope ∂MOC/∂*τ*_{0} gives a sense of how strongly each cell depends on the wind. These values are given in the first column of Table 2, along with the value of *R*^{2} for the regression. The *R*^{2} values reveal that the linear fit is very good in most cases.

#### 1) Fixed flux boundary condition

The MOC transports are rather insensitive to the wind in the fixed flux experiments. Here, MOC_{lower} shows no correlation with *τ*_{0}, varying in a narrow range about 0.4 Sv; MOC_{upper} is quite weak for the weakest winds (*τ*_{0} = 0.05 and 0.1 N m^{−2}), but for the rest of the experiments (0.125 N m^{−2} ≥ *τ*_{0} ≥ 0.3 N m^{−2}), the changes in MOC_{upper} are slight: it increases only from 0.5 to 0.6 Sv over this range. The linear fit for ∂MOC_{upper}/∂*τ*_{0} (Table 2) shows a sensitivity ¼ of that of the Ekman circulation. Examination of the structure of Ψ_{iso} show that, for weak winds, the upper cell becomes confined more and more to the surface diabatic layer and does not reach the interior. Because MOC_{lower} does not change, MOC_{upwell} follows the changes in MOC_{upper}. The small changes in residual MOC reflect the fact that, as the magnitude of increases with the wind, Ψ_{eddy} also strengthens (becoming increasingly negative), leading to a high degree of compensation between mean and eddy circulations.

#### 2) Relaxation boundary condition

We implement the interactive boundary condition in the MITgcm by relaxing the temperature in the top model level, referred to as *T _{s}*, to a prescribed function of latitude

*T**(

*y*). For the base-case winds (

*τ*

_{0}= 0.2 N m

^{−2}), we wish to have the same effective heat flux as the fixed-flux reference case described above. For a layer of depth Δ

*z*subject to relaxation at a rate

*λ*, the effective heat flux is

We chose a relaxation time scale of *λ* = 30 day^{−1}.^{1} Given *T _{s}* from the base-case fixed-flux experiment, the desired heat flux

*Q*[Eq. (1)], and Δ

*z*= 10 m, this expression can be rearranged to find

*T**. As expected, when the

*τ*

_{0}= 0.2 N m

^{−2}simulation is run with this forcing, it reaches the same equilibrium as the base-case fixed-flux state described in the previous section, with the same MOC transport, because the heat flux felt by the ocean is nearly unchanged. However, when

*τ*

_{0}is changed,

*T*can and does change, resulting in an altered air–sea heat flux and, evidently, greater sensitivity of Ψ

_{s}_{iso}to the winds.

The results of these experiments are also shown in Fig. 5. The changes are significantly larger than the fixed-flux case. Both MOC_{upper} and MOC_{lower} increase with stronger winds; this means a strengthening of the upper cell (because it is positive: i.e., clockwise) and a weakening of the lower (negative, counterclockwise) cell. The linear fit (Table 2) shows that MOC_{upper} is nearly twice as sensitive as the fixed-flux case. Because the changes in are the same under both boundary conditions, the higher sensitivity implies that the magnitude Ψ_{eddy} is not as sensitive to *τ*_{0}, leading to less compensation.

Viebahn and Eden (2010) performed a very similar experiment, simulating only an upper cell and using a relaxation boundary condition for buoyancy. Their results are broadly consistent with ours: a sensitivity of the residual circulation much weaker than the sensitivity of the Ekman circulation. However, they observed decreasing sensitivity with increasing winds, whereas the trend in our MOC_{upper} appears quite linear. This qualitative difference is most likely attributable to the different northern boundary; they had a small, unforced basin attached to the northern edge of their channel, rather than a sponge layer.

## 3. The surface buoyancy boundary condition

Our experiments make it clear that a residual overturning driven by a fixed buoyancy flux is less sensitive to the winds than one with an interactive buoyancy flux. In this section, we seek to understand this behavior diagnostically through the residual buoyancy budget, using the framework of MR03.

### a. Transformed Eulerian-mean buoyancy budget

We begin by reviewing some essential elements of TEM theory (Andrews and McIntyre 1976; Andrews et al. 1987; Treguier et al. 1997; Plumb and Ferrari 2005). The reader is referred to MR03 for a complete discussion of the theory in the context of ACC dynamics.

The time and zonally averaged buoyancy equation for our domain (outside of the sponge layer) is

where is the downward buoyancy flux from the surface forcing (we have neglected the relatively small fluxes due to diffusion). The goal of TEM theory is to simplify the eddy-flux terms by separating them into advective and diabatic components. The eddy advection can then be combined with the mean advection in a residual streamfunction Ψ_{res},

MR03 choose to define the TEM eddy streamfunction as

The TEM residual streamfunction Ψ_{res}, defined using (11) and (12), is nearly identical to Ψ_{iso}. Now, (10) can be manipulated into the form

where the Jacobian term *J* represents advection by the residual circulation and the factor *μ* measures the diabatic eddy flux contribution,

When the eddy flux is directed along mean isopycnals, *μ* = 1 and the second term on the RHS of (13) vanishes. If both terms on the RHS are zero, as expected in the ocean interior, then means that Ψ_{res} is constant along isopycnals.

At this point, MR03 make several assumptions to arrive at an analytic solution. First, they assume the existence of a mixed layer of fixed depth *h _{m}* in which . Following Treguier et al. (1997), they assume that

*μ*= 1 in the ocean interior and varies from 1 at the base of the mixed layer to 0 at the surface. The buoyancy flux

*B*is also assumed to reach zero by the base of the mixed layer. Integrating (13) over this mixed layer, one obtains

where is the surface mixed layer buoyancy. We have also defined

as the mixed layer–integrated diabatic eddy flux divergence. Equation (15) states that advection by the residual flow across the mixed layer buoyancy gradient is balanced by diabatic forcing and diabatic eddy fluxes.

In the MR03 model [of which (15) is a central component], Ψ_{res} reaches its full value at the base of the mixed layer. However, in our model, the surface diabatic layer (200–300 m) extends much deeper than the shallow mixed layer (~50 m). Figure 3 makes it clear that Ψ_{iso} (approximately equivalent to Ψ_{res}) does not reach its full interior value until below the diabatic layer. To overcome this complication, we express *B*, *D*, and as functions of the vertically averaged buoyancy within the diabatic layer *b*_{dl}. We can then write (15) approximately as

We can see by plotting the three terms of this equation (Fig. 6) that the agreement is good throughout most of the domain. This diagnosis shows that the diabatic forcing *B* determines the strength and sense of the interior MOC. The diabatic eddy-flux term *D* is small but not negligible; it generally opposes *B*, resulting in a weaker Ψ_{iso}. The largest imbalance in this approximate form arises in the region associated with uppermost counterclockwise cell (between *T* = 5°C and *T* = 6°C), which contains a recirculation cell entirely within the diabatic layer, a complication not considered in the MR03 theory. This surface cell is not the focus of our analysis.

### b. Buoyancy-flux sensitivity to winds

The residual buoyancy budget as expressed by (15) or (17) already reveals the strong constraint imposed on the MOC by a fixed surface buoyancy flux: because the term *B* cannot change, changes in the MOC must be accompanied by changes in or *D*. In the relaxation case, in contrast, *B* can also change, implying a higher degree of freedom for the MOC. This freedom is reflected in the higher MOC sensitivity in the relaxation experiments.

As described above for the reference case, we can diagnose the forcing terms *B* and *D* from each of our experiments to understand how these terms change with the wind; this is shown in Fig. 7, which contains contour plots of *B* and *D* as functions of *y* and *τ*_{0}. Also plotted are contours of the zonal-mean SST, from which it is easy to see the changes in . From this figure, we can see that two factors contribute to the strengthening of the upper cell in the fixed-flux case. First, the diabatic eddy flux *D* (Fig. 7c), which generally opposes the heating from *B* centered on *y* = 666 km, decreases with increasing winds, leading to greater total buoyancy gain in this region. Second, the SST contours in this region spread apart as the winds increase, decreasing , which further contributes in the increase in Ψ_{res}. On the other hand, in the lower-cell formation region (the most southern part of the domain), trends in *D* and evidently cancel, leading to no trend in the Ψ_{res} associated with the lower cell.

The changing air–sea buoyancy flux *B* in the relaxation case is evident in Fig. 7b. The flux is everywhere increasing as the winds increase, in accord with the fact that SSTs are decreasing [see (9); SSTs also change in the fixed-flux case; however, because the flux is not interactive, this has no effect on *B*]. This is completely consistent with the increased upper-cell transport and decreased lower-cell transport. In comparison with the fixed-flux case, the changes in are less significant, resulting from the fact that SST is being relaxed to the same function of *y* in all experiments. Changes in *D* seem insignificant for the upper-cell region but still potentially important for the lower cell.

Dependence of the air–sea buoyancy flux *B* on wind stress was observed by Badin and Williams (2010) in a similar yet coarse-resolution model. Their study also noted the sensitivity of *B* to the choice of Gent–McWilliams eddy-transfer coefficient. In our interactive buoyancy-flux experiments, both the eddy transfer and the buoyancy flux are free to respond to changing winds, resulting in a tangled equilibration problem. The diagnostics presented in this section merely show how the buoyancy budget is consistent with the residual circulation; they do not explain the magnitude of the sensitivity. For that, we need to look closer at the eddy circulation itself.

## 4. Constraints on the eddy circulation

In this section, we seek to understand what sets the strength of the eddy circulation. This discussion is most relevant to the interactive buoyancy-flux experiments, whose residual circulation cannot be assumed a priori based on knowledge of the buoyancy flux. The essential question is, how well can we estimate the sensitivities of Ψ_{iso} reported in Table 2 based on first principles?

### a. Decomposing the eddy circulation: Slope and diffusivity

In the adiabatic interior, the TEM eddy circulation can be written as

This form is identical to the earlier definition (12) when *μ* = 1, which is a good approximation away from the surface diabatic layer and sponge layer and very close to Ψ_{eddy} as defined in (8).

Assuming a flux gradient relationship , where *K* is the eddy diffusivity, we can write (18) as

where is the mean isopycnal slope. Equation (19) is the basis of the famous Gent–McWilliams parameterization for mesoscale eddies (Gent and McWilliams 1990; Gent et al. 1995). Here, it is simply a rearrangement of the definition of Ψ* given the definition of *K*. Using the definition of [Eq. (6)], the residual circulation then becomes

This expression is a centerpiece of the MR03 model.

Viebahn and Eden (2010) applied (20) to their eddy-resolving model in order to ascertain the relative importance of changes in *K* and *s*. We follow a very similar path. Consider a reference state in which : the variables for this state will be denoted *K*_{ref}, *s*_{ref}, etc. For different values of *τ*_{0}, the departures of these variables from the reference state will be expressed as Δ*K*, where Δ*K* = *K* − *K*_{ref}, and similarly for the other variables. Using this notation, we can express Ψ_{res} for any *τ*_{0} state as

where a quadratic Δ term has been dropped (note that and ). The first term in (22) expresses the linear scaling of the Ekman-driven circulation with the wind. The second term expresses the eddy response. If we can develop a theory for the fractional changes in *K* and *s* with changing *τ*_{0}, we can effectively predict the MOC departure from a reference state for a change in winds.

Viebahn and Eden (2010) found that changes in *s* were very small compared to changes in *K* and that the changes in Ψ* could therefore be attributed primarily to changes in *K*. To test this idea in our model, we calculate , Δ*K*/*K*_{ref}, and Δ*s*/*s*_{ref} from the model output. The calculation is performed at a depth of 477 m, below the surface diabatic layer but shallow enough to see all the MOC cells [above this depth, we find that (18) is not a very good approximation of Ψ_{iso}]. The terms are plotted in Fig. 8 as a function of *y* and *τ*_{0}. We see that Ψ* changes by about 50% from the reference case in either direction (weaker or stronger winds), and *K* undergoes changes in magnitude almost as large. However, *s* is notably less sensitive, weakening by 20% for weak winds and barely changing at all for stronger winds. Changes in *s* are most significant in the southernmost part of the domain, where new isopycnals outcrop with increasing winds. For the fixed-flux experiments, in contrast, *s* undergoes large changes in a wider part of the domain (not shown).

The relative insensitivity of *s* seems somewhat inevitable given the boundary conditions. Because the buoyancy is relaxed to prescribed values at both the surface and the northern boundary, the large-scale isopycnal slope is effectively prescribed as well (of course, small changes in surface buoyancy are necessary to bring about changes in heat flux, as seen in Fig. 7). Only isopycnals that do not outcrop are unconstrained on the southern edge, resulting in higher values of Δ*s*/*s*_{ref} in the far southern part of the domain. Viebahn and Eden (2010) found Δ*s* to be small in an experiment with no sponge layer.

Focusing on the same surface (*z* = −477 m), we can use (22) to directly estimate MOC_{upper} and MOC_{lower} by picking points in *y* that correspond with the maximum and minimum values of Ψ_{iso} (these points do not move significantly in space with changes in *τ*_{0}). By calculating Δ*s* and Δ*K* at these points, we can evaluate (22). The linear MOC sensitivities produced in this way are given in the third column of Table 2. These sensitivities agree very well with the values given by Ψ_{iso}, indicating that (22) is a good approximation.

Given the observed smallness of Δ*s*, we can ask, to what extent is the sensitivity of the MOC due to Δ*K*? To answer this question, we evaluate (22) with Δ*s* = 0 and compute the linear sensitivity. For comparison, we also do the opposite, setting Δ*K* = 0 and using only Δ*s*. The results, given in the fourth column of Table 2, indicate that Δ*K* is the dominant factor in the upper-cell sensitivity in both the fixed-flux and relaxation experiments. Especially in the relaxation experiment, the sensitivity due to Δ*s* alone is close to the sensitivity, suggesting a negligible role for Δ*s*. In contrast, Δ*K* and Δ*s* seem to play equal roles in the lower-cell sensitivity.

### b. Eddy diffusivity dependence on wind stress

Given the prominent role of Δ*K* in determining the MOC sensitivity, we focus now on understanding its scaling behavior with the winds. As a starting point, we plot the full *K*(*y*, *z*) for three different values of *τ*_{0} in Fig. 9. In general, *K* is positive nearly everywhere and appears intensified very near the surface and toward the bottom, with a minimum at middepth. The details of the vertical structure of *K* are interesting but are not our focus here (a paper on this topic is in preparation). For now, we simply note that the spatial structure does not change qualitatively with *τ*_{0}, allowing us to imagine a fixed spatial structure that simply scales with *τ*_{0} (Viebahn and Eden 2010 found a strikingly similar spatial pattern). Many studies, including MR03 and Visbeck et al. (1997), have assumed that *K* itself is proportional to *s*. Instead, we employ a mixing-length theory, which relates K to the eddy kinetic energy and thus to the mechanical energy balance.

Mixing-length theory (Taylor 1921; Prandtl 1925) claims that the eddy diffusivity can be expressed as a characteristic eddy velocity *V _{e}* times an eddy length scale

*L*, such that . Many authors have applied this idea to estimate eddy diffusivities in the ocean (Holloway 1986; Keffer and Holloway 1988; Davis 1991; Stammer 1998; Eden and Greatbatch 2008; Ferrari and Nikurashin 2010).

_{e}^{2}A general theory predicting

*V*and

_{e}*L*for geostrophic turbulence does not yet exist, but the topic is a very active area of research (Held and Larichev 1996; Lapeyre and Held 2003; Thompson and Young 2007).

_{e}Cessi (2008) suggested that the appropriate *V _{e}* to use for the buoyancy diffusivity is the barotropic eddy velocity (i.e., the RMS anomaly of the vertically averaged velocity), because barotropic stirring can most efficiently mix buoyancy across sloping isopycnals. We make the key assumption that this value and thus

*K*itself are proportional (but not necessarily equal) to the bottom eddy velocity. In terms of eddy kinetic energy (the square of

*V*), this statement becomes

_{e}where is the velocity anomaly at the bottom. The angle brackets indicate an average in *x*, *y* and time. (The ideas could be easily extended to include dependence on *y*, but here we find it simpler to concentrate on the domain average.)

Following Cessi et al. (2006) and Cessi (2008), we consider the mechanical energy budget. Because our model employs linear bottom drag, the leading-order balance of the system is

with additional small contributions from viscous dissipation, side drag, and conversion to potential energy (the zonal-mean surface velocity has been used because *τ _{s}* is constant in

*x*and time). Physically, (24) expresses the fact that the wind power input to the system is dissipated by bottom drag. Cessi et al. (2006) found this balance to hold well in a similar numerical model. We checked (24) in our model and found it to hold not only globally but also in a zonal-average budget to within 10% error (not shown).

The bottom velocity **u*** _{b}* includes both the mean and eddy velocity, . The term is negligible in comparison to . Furthermore, we already know from (5) that . This allows (24) to be rearranged to the form

It is important to note that the eddy energy in this expression depends only on the baroclinic shear . The large barotropic velocity due to the absence of topography does not affect the eddy energy balance. Furthermore, because topographic form drag does not participate in the energy cycle (Ferrari and Wunsch 2009), we can expect (25) to hold in more realistic models with topography. (Note that in the presence of topopgraphy, .) The baroclinic shear can be obtained from the thermal-wind equation,

Because the large-scale meridional buoyancy gradient is determined by the large-scale forcing, we can expect this thermal-wind contribution to remain approximately constant.^{3} This suggests the scaling relationship

where Δ*T* is the large-scale temperature difference across the channel, set by the relaxation SST. We have also used the fact that 〈*τ _{s}*〉 = 2

*τ*

_{0}/π. Our mixing-length hypothesis claims that

*K*is related to this quantity as , where the mixing-length constant

*L*has absorbed the unknown constant of proportionality between EKE

_{e}_{bt}and EKE

*.*

_{b}We can diagnose all these quantities from the model to test our ideas. We assume that the *L _{e}* = constant = 30 km. In reality, the mixing length also varies by ~10% (as diagnosed from the simulations), but we can achieve decent agreement without considering these effects, and a theory for

*L*is beyond the scope of this paper. The quantities EKE

_{e}*, EKE*

_{b}_{bt}, (〈

*K*〉/

*L*)

_{e}^{2}, and the scaling prediction from (27) are all plotted in Fig. 10 on a logarithmic scale as a function of

*τ*

_{0}. Both the fixed-flux and relaxation cases are plotted. The three diagnosed quantities and the theoretical prediction show similar slopes, especially for high value of

*τ*

_{0}. The small departures of 〈

*K*〉 from the EKE values can be explained by our neglect of changes in

*L*. The small departures of EKE

_{e}*from the scaling theory likewise can be explained by our neglect of second-order energy sources and sinks in (24). However, based on the general agreement, we conclude that a useful approximation for the eddy buoyancy diffusivity in our model is*

_{b}When applying this formula locally in space, we must, however, expect errors due to the changing spatial structure of *K* shown in Fig. 9.

### c. Predicting the MOC sensitivity

Using the scaling from (28) in (22), along with the assumption that changes in *s* can be neglected to first order, we arrive at

From this equation, the linear sensitivities of MOC_{upper} and MOC_{lower} can be calculated analytically given and which are given in the final column of Table 2. Two important points must be kept in mind in interpreting these values. First, such estimates can only be as good as the Δ*K*-only sensitivity already presented, which we noted was most accurate for the relaxation-case upper cell. Second, because and are nearly the same for both fixed-flux and relaxation experiments, the scaling for Δ*K* produces nearly identical predictions for the different cases. The final prediction from the scaling theory of the sensitivity of the upper cell in the relaxation case, 6.9 Sv (N m^{−2})^{−1}, is higher than the Δ*K*-only sensitivity, 5.9 Sv (N m^{−2})^{−1}, which itself is higher than the true sensitivity, 4.5 Sv (N m^{−2})^{−1}. However, all these values are significantly weaker than the sensitivity of , 11.1 Sv (N m^{−2})^{−1}. The agreement between the relaxation-case lower-cell sensitivity from the scaling and the true sensitivity is spurious, a case of two wrongs [the neglect of Δ*s* and the failure of (28) locally for the lower-cell region] making a right. Although our scaling theory is far from comprehensive, we are encouraged by the agreement for the upper-cell relaxation case and consider it a useful stepping stone in a difficult problem.

## 5. Discussion and conclusions

One important conclusion of this study is that the sensitivity of the Southern Ocean MOC to the winds depends on the surface boundary condition for buoyancy. This is not an immediately intuitive result, because the winds are a purely mechanical forcing. However, it becomes clear once one considers the TEM (or, equivalently, isopycnal average) point of view expressed by (15): in a quasi-adiabatic ocean interior, the residual MOC is primarily set by diabatic water-mass transformation at the surface, and, if the winds are unable to alter the transformation rates (as in the fixed-buoyancy-flux case), the sensitivity of the MOC is weak. In fact, evidence of this point emerges from the existing literature when comparing different models. For instance, Hallberg and Gnanadesikan (2006) used a predominantly fixed-flux surface boundary condition and found a relatively weak sensitivity of the residual MOC to increased winds. In contrast, Wolfe and Cessi (2010) used a relaxation boundary condition and found much greater sensitivity; in certain locations, they found an increase in residual MOC transport almost equal to the increase in Ekman transport, the upper limit of the sensitivity. The increased transport was accompanied by increased transformation in both southern and northern high latitudes. Although our model contains only an ACC channel, it manages to qualitatively reproduce the behavior of both these two different models just by changing the surface boundary condition. Similar conclusions were reached by Bugnion et al. (2006), using an adjoint method in a coarse-resolution model, and by Badin and Williams (2010).

The surface boundary condition of the real ocean is mixed. Certain contributions to the air–sea buoyancy flux, such as net shortwave radiation and precipitation, are largely independent of the SST and surface winds. Latent and sensible heat fluxes, on the other hand, are interactive (Haney 1971). For the winds to play a strong role in modulating the residual MOC, as envisioned by Toggweiler and Russell (2008), our study suggests that the interactive fluxes must dominate. It should therefore be a top priority to continue to improve our understanding of the processes that determine the air–sea buoyancy flux in the Southern Ocean—including sea ice processes, which we have completely neglected—and whether these components are sensitive to changes in wind or other climate changes.

Of the various simplifications we have made, perhaps the most restrictive and unrealistic is the fixed stratification imposed by the northern boundary sponge layer. In fact, many of the related studies we have cited have focused explicitly on the question of what sets the stratification (MR03; Henning and Vallis 2005; Wolfe and Cessi 2010). In the analytical model of MR03, the thermocline depth was found to be proportional to , but their eddy closure (*K* ∝ |*s*|) does not hold in our eddy-resolving model. Henning and Vallis (2005) found a weaker scaling of the stratification with the wind in an eddy-resolving model of a channel coupled to a basin. Such results are encouraging because, if the stratification dependence on *τ*_{0} is weak, it is more reasonable to approximate it as fixed, as we have done. Nevertheless, tests of our results in more realistic global eddy-resolving models are required.

Finally, we developed a scaling theory for the eddy diffusivity and used it estimate the MOC sensitivity. Traditionally, scaling theories for eddies have been based on ideas from linear baroclinic instability, and the eddy diffusivity is assumed to be somehow proportional to the isopycnal slope (Green 1970; Stone 1972; Killworth 1997; Visbeck et al. 1997). Although baroclinic instability plays a crucial role in the energy cycle of our model, linear theory cannot predict the fully equilibrated eddy energy. Instead we have followed some of the ideas developed by Cessi (2008), invoking the mechanical energy balance to gain insight into the eddy energy and diffusivity. Consequently, our scaling theory for the eddy diffusivity (27) includes a dependence on both the wind stress parameter *τ*_{0} and the bottom drag *r _{b}* but not on the isopycnal slope

*s*. The scaling shows good agreement with the GCM results. Furthermore, we think it represents a promising way forward in understanding the role of eddies in the equilibration of the Southern Ocean.

We have examined only steady states, but the time-dependent response to wind changes is important and interesting. Meredith and Hogg (2006) have suggested Southern Ocean eddies can respond very fast (~1 yr) to changes in wind, whereas Treguier et al. (2010) found that the interannual MOC variability in a realistic model was dominated by Ekman transport, with little eddy compensation. This issue deserves further study as well.

## Acknowledgments

Raf Ferrari and Malte Jansen made several helpful suggestions. Comments from Paola Cessi, two anonymous reviewers, and the editor greatly improved the manuscript. The observational wind and buoyancy flux data are from the Research Data Archive (RDA), which is maintained by the Computational and Information Systems Laboratory (CISL) at the National Center for Atmospheric Research (NCAR). NCAR is sponsored by the National Science Foundation (NSF). The original data are available from the RDA (http://dss.ucar.edu) in dataset ds260.2.

## REFERENCES

_{2}

_{2}change

## Footnotes

^{1}

This choice of parameters corresponds to a sensitivity of ∂*Q*_{eff}/∂*T _{s}* ~ 15 W m

^{−2}K

^{−1}.

^{2}

Ferrari and Nikurashin (2010) recently refined the idea to include the modulation of *L _{e}* by the presence of mean flows, and there is indeed mounting evidence that the spatial variations in

*K*in the Southern Ocean are modulated by the strong jets found there (Marshall et al. 2006; Smith and Marshall 2009; Abernathey et al. 2010; Naveira Garabato et al. 2011).

^{3}

This is equivalent to assuming that the baroclinic transport in the model is “saturated” (Straub 1993), which is indeed the case for our experiments.