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 CO2 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.

Fig. 1.

Maps of the observed surface forcing in the Southern Ocean, averaged from the Common Ocean Reference Experiment (CORE2) dataset over the period 1949–2006 (Large and Yeager 2009): (left) The wind stress in N m−2, with the magnitude indicated by the colored shading and the direction by the arrows and (right) the buoyancy-equivalent heat flux in W m−2, which includes contributions from longwave and shortwave radiative fluxes; latent and sensible heat fluxes; and the buoyancy fluxes due to evaporation, precipitation, and runoff.

Fig. 1.

Maps of the observed surface forcing in the Southern Ocean, averaged from the Common Ocean Reference Experiment (CORE2) dataset over the period 1949–2006 (Large and Yeager 2009): (left) The wind stress in N m−2, with the magnitude indicated by the colored shading and the direction by the arrows and (right) the buoyancy-equivalent heat flux in W m−2, which includes contributions from longwave and shortwave radiative fluxes; latent and sensible heat fluxes; and the buoyancy fluxes due to evaporation, precipitation, and runoff.

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.

Table 1.

Parameters used in the numerical model reference experiment.

Parameters used in the numerical model reference experiment.
Parameters used in the numerical model reference experiment.

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

 
formula

and Q = 0 north of this point, with Q0 = 10 W m−2. The term Ly 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

 
formula

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 = Lsponge) to 7 day−1 at the northern boundary (y = Ly). 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

 
formula

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

 
formula

where rb is a bottom drag coefficient and ub the horizontal component of the bottom velocity.

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.

Fig. 2.

A 3D snapshot of the model’s temperature field, revealing the mesoscale eddy field. The temperatures range from 0° to 8°C. Overlaid on top are depictions of the wind stress and heat-flux surface forcing. To the right is the zonal- and time-mean zonal velocity , which ranges from 0 to 25 cm s−1. The contour interval for is 2.5 cm s−1. Overlaid in white are the 1°, 3°, and 5°C isotherms.

Fig. 2.

A 3D snapshot of the model’s temperature field, revealing the mesoscale eddy field. The temperatures range from 0° to 8°C. Overlaid on top are depictions of the wind stress and heat-flux surface forcing. To the right is the zonal- and time-mean zonal velocity , which ranges from 0 to 25 cm s−1. The contour interval for is 2.5 cm s−1. Overlaid in white are the 1°, 3°, and 5°C isotherms.

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

 
formula

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 ≡ 106 m3 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

 
formula

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

 
formula

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.

Fig. 3.

The residual MOC streamfunction Ψiso (left) as originally diagnosed in isopycnal coordinates and (right) mapped back to depth coordinates. The units are Sv, and the contour interval is 0.1 Sv. The solid black line in (left) indicates the mean SST, and the gray lines are the 5% [in (left) only] and 95% levels of the SST CDF. The dotted black is the southern boundary of the sponge layer. The mean T contours are also shown in (right) in black, and the contour interval is 0.5°C.

Fig. 3.

The residual MOC streamfunction Ψiso (left) as originally diagnosed in isopycnal coordinates and (right) mapped back to depth coordinates. The units are Sv, and the contour interval is 0.1 Sv. The solid black line in (left) indicates the mean SST, and the gray lines are the 5% [in (left) only] and 95% levels of the SST CDF. The dotted black is the southern boundary of the sponge layer. The mean T contours are also shown in (right) in black, and the contour interval is 0.5°C.

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,

 
formula

Both and Ψeddy are shown in Fig. 4. Their magnitudes are large, but they oppose each other, leaving Ψiso as a small residual. Because the dependence of on the wind is clear from (6), the difficulty in understanding the residual MOC sensitivity to the winds lies in Ψeddy.

Fig. 4.

(left) The Eulerian-mean streamfunction and (right) the eddy streamfunction Ψeddy, as defined by (8). The units are Sv, and the contour interval is 0.5 Sv. Otherwise, as in Fig. 3 (right).

Fig. 4.

(left) The Eulerian-mean streamfunction and (right) the eddy streamfunction Ψeddy, as defined by (8). The units are Sv, and the contour interval is 0.5 Sv. Otherwise, as in Fig. 3 (right).

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 MOCupper and MOClower. Besides the individual upper and lower cells, there is a third relevant quantity: the total volume flux of upwelled deep water, MOCupwell = MOCupper − MOClower. This value is also shown in Fig. 5, along with the strength of , the Ekman circulation. In almost all cases, MOCupper, MOClower, and MOCupwell are weaker than .

Fig. 5.

A summary of the MOC cell strength in all of the different experiments. The Ekman circulation is shown in black, and the residual circulations of the various MOC cells (upper, lower, and net upwelling) are plotted in color. Fixed-surface-flux experiments are represented in blue; surface-relaxation experiments are represented in orange. The shapes correspond to the values of MOClower, MOCupper, and MOCupwell. The reference case, τ0 = 0.2 N m−2, is indicated by the dotted line.

Fig. 5.

A summary of the MOC cell strength in all of the different experiments. The Ekman circulation is shown in black, and the residual circulations of the various MOC cells (upper, lower, and net upwelling) are plotted in color. Fixed-surface-flux experiments are represented in blue; surface-relaxation experiments are represented in orange. The shapes correspond to the values of MOClower, MOCupper, and MOCupwell. The reference case, τ0 = 0.2 N m−2, is indicated by the dotted line.

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 R2 for the regression. The R2 values reveal that the linear fit is very good in most cases.

Table 2.

Linear MOC dependence on wind ∂MOC/∂τ0, as determined by least squares fit. The value of R2 for the linear regression is given in parenthesis, a measure of the goodness of fit. The values are computed at fixed points in space near where maxima and minima of Ψiso occur: z = −477 m, y = 1150 km (upper cell), and y = 300 km (lower cell). The first column shows Ψiso, and the second column shows . The rightmost four columns represent the approximations produced by (22), (22) with Δs set to zero, (22) with ΔK set to zero, and finally (29), the prediction for the MOC sensitivity given by neglecting Δs and and assuming that K scales locally with (28).

Linear MOC dependence on wind ∂MOC/∂τ0, as determined by least squares fit. The value of R2 for the linear regression is given in parenthesis, a measure of the goodness of fit. The values are computed at fixed points in space near where maxima and minima of Ψiso occur: z = −477 m, y = 1150 km (upper cell), and y = 300 km (lower cell). The first column shows Ψiso, and the second column shows . The rightmost four columns represent the approximations produced by (22), (22) with Δs set to zero, (22) with ΔK set to zero, and finally (29), the prediction for the MOC sensitivity given by neglecting Δs and and assuming that K scales locally with (28).
Linear MOC dependence on wind ∂MOC/∂τ0, as determined by least squares fit. The value of R2 for the linear regression is given in parenthesis, a measure of the goodness of fit. The values are computed at fixed points in space near where maxima and minima of Ψiso occur: z = −477 m, y = 1150 km (upper cell), and y = 300 km (lower cell). The first column shows Ψiso, and the second column shows . The rightmost four columns represent the approximations produced by (22), (22) with Δs set to zero, (22) with ΔK set to zero, and finally (29), the prediction for the MOC sensitivity given by neglecting Δs and and assuming that K scales locally with (28).

1) Fixed flux boundary condition

The MOC transports are rather insensitive to the wind in the fixed flux experiments. Here, MOClower shows no correlation with τ0, varying in a narrow range about 0.4 Sv; MOCupper 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 MOCupper are slight: it increases only from 0.5 to 0.6 Sv over this range. The linear fit for ∂MOCupper/∂τ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 MOClower does not change, MOCupwell follows the changes in MOCupper. 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 Ts, 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

 
formula

We chose a relaxation time scale of λ = 30 day−1.1 Given Ts 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, Ts can and does change, resulting in an altered air–sea heat flux and, evidently, greater sensitivity of Ψ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 MOCupper and MOClower 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 MOCupper 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 MOCupper 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

 
formula

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,

 
formula

MR03 choose to define the TEM eddy streamfunction as

 
formula

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

 
formula

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

 
formula

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 hm 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

 
formula

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

 
formula

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 bdl. We can then write (15) approximately as

 
formula

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.

Fig. 6.

The terms in the approximate form of the Marshall–Radko balance (17). The plot is shown as a function of T on the bottom of the x axis, but it can also be considered a function of y, whose corresponding values are shown at the top of the x axis.

Fig. 6.

The terms in the approximate form of the Marshall–Radko balance (17). The plot is shown as a function of T on the bottom of the x axis, but it can also be considered a function of y, whose corresponding values are shown at the top of the x axis.

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.

Fig. 7.

The forcing terms of the surface residual buoyancy budget (22) for changing values of wind τ0, expressed in units of W m−2 equivalent by multiplying by ρ0cp()−1. Shown are (a) the air–sea buoyancy flux B for the fixed-flux case and (c) the diabatic eddy flux D. (b),(d) As in (a),(c), but for the relaxation surface boundary condition. The thin black lines are contours of the zonal-average SST, and the contour interval is 0.5°C, from which changes in the surface buoyancy gradient ∂bs/∂y can be inferred. The thick dashed black lines indicate the boundaries of the regions of applied surface heating and cooling from the reference experiment τ0 = 0.2 N m−2.

Fig. 7.

The forcing terms of the surface residual buoyancy budget (22) for changing values of wind τ0, expressed in units of W m−2 equivalent by multiplying by ρ0cp()−1. Shown are (a) the air–sea buoyancy flux B for the fixed-flux case and (c) the diabatic eddy flux D. (b),(d) As in (a),(c), but for the relaxation surface boundary condition. The thin black lines are contours of the zonal-average SST, and the contour interval is 0.5°C, from which changes in the surface buoyancy gradient ∂bs/∂y can be inferred. The thick dashed black lines indicate the boundaries of the regions of applied surface heating and cooling from the reference experiment τ0 = 0.2 N m−2.

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

 
formula

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

 
formula

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

 
formula

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 Kref, sref, etc. For different values of τ0, the departures of these variables from the reference state will be expressed as ΔK, where ΔK = KKref, and similarly for the other variables. Using this notation, we can express Ψres for any τ0 state as

 
formula
 
formula

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/Kref, and Δs/sref 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).

Fig. 8.

Fractional changes in (a) eddy circulation (b) eddy diffusivity ΔK/Kref, and (c) isopycnal slope Δs/sref from the reference case (indicated by the dashed black line) computed at 477-m depth. The black contours are the mean isotherms at this depth, and the contour interval is 0.5°C.

Fig. 8.

Fractional changes in (a) eddy circulation (b) eddy diffusivity ΔK/Kref, and (c) isopycnal slope Δs/sref from the reference case (indicated by the dashed black line) computed at 477-m depth. The black contours are the mean isotherms at this depth, and the contour interval is 0.5°C.

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/sref 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 MOCupper and MOClower 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.

Fig. 9.

Flux-gradient buoyancy diffusivity K(y, z) for (left to right) three different wind strengths and fixed-flux boundary condition, shown with contour interval 500 m2 s−1. The black contours are the mean isotherms, and the contour interval is 0.5°C.

Fig. 9.

Flux-gradient buoyancy diffusivity K(y, z) for (left to right) three different wind strengths and fixed-flux boundary condition, shown with contour interval 500 m2 s−1. The black contours are the mean isotherms, and the contour interval is 0.5°C.

Mixing-length theory (Taylor 1921; Prandtl 1925) claims that the eddy diffusivity can be expressed as a characteristic eddy velocity Ve times an eddy length scale Le, 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).2 A general theory predicting Ve and Le 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).

Cessi (2008) suggested that the appropriate Ve 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 Ve), this statement becomes

 
formula

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

 
formula

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 ub 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

 
formula

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,

 
formula

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

 
formula

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 Le has absorbed the unknown constant of proportionality between EKEbt and EKEb.

We can diagnose all these quantities from the model to test our ideas. We assume that the Le = 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 Le is beyond the scope of this paper. The quantities EKEb, EKEbt, (〈K〉/Le)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 Le. The small departures of EKEb 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

 
formula

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

Fig. 10.

Globally averaged EKEs diagnosed from the model. The quantities plotted are the barotropic EKE (square), the bottom EKE (circle), and the EKE implied by the diffusivity K (triangle), assuming a constant mixing length of 30 km. The black line is the EKE predicted by the scaling relation (27). The fixed-heat-flux experiments are white, whereas the relaxation experiments are gray.

Fig. 10.

Globally averaged EKEs diagnosed from the model. The quantities plotted are the barotropic EKE (square), the bottom EKE (circle), and the EKE implied by the diffusivity K (triangle), assuming a constant mixing length of 30 km. The black line is the EKE predicted by the scaling relation (27). The fixed-heat-flux experiments are white, whereas the relaxation experiments are gray.

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

 
formula

From this equation, the linear sensitivities of MOCupper and MOClower 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 rb 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

REFERENCES
Abernathey
,
R.
,
J.
Marshall
,
E.
Shuckburgh
, and
M.
Mazloff
,
2010
:
Enhancement of mesoscale eddy stirring at steering levels in the Southern Ocean
.
J. Phys. Oceanogr.
,
40
,
170
185
.
Anderson
,
R. F.
,
S.
Ali
,
L. I.
Bradtmiller
,
S. H. H.
Nielsen
,
M. Q.
Fleisher
,
B. E.
Anderson
, and
L. H.
Burckle
,
2009
:
Wind-driven upwelling in the Southern Ocean and the deglacial rise in atmospheric CO2
.
Science
,
323
,
1143
1150
.
Andrews
,
D. G.
, and
M. E.
McIntyre
,
1976
:
Planetary waves in horizontal and vertical shear: The generalized Eliassen-Palm relation and the mean zonal acceleration
.
J. Atmos. Sci.
,
33
,
2031
2058
.
Andrews
,
D. G.
,
J. R.
Holton
, and
C. B.
Leovy
,
1987
:
Middle Atmosphere Dynamics
.
International Geophysical Series, Vol. 40, Academic Press, 489 pp
.
Badin
,
G.
, and
R. G.
Williams
,
2010
:
On the buoyancy forcing and residual circulation in the Southern Ocean: The feedback from Ekman and eddy transfer
.
J. Phys. Oceanogr.
,
40
,
295
311
.
Bugnion
,
V.
,
C.
Hill
, and
P. H.
Stone
,
2006
:
An adjoint analysis of the meridional overturning circulation in an ocean model
.
J. Climate
,
19
,
3732
3751
.
Cenedese
,
C.
,
J.
Marshall
, and
J. A.
Whitehead
,
2004
:
A laboratory model of thermocline depth and exchange fluxes across circumpolar fronts
.
J. Phys. Oceanogr.
,
34
,
656
668
.
Cerovečki
,
I.
,
R. A.
Plumb
, and
W.
Heres
,
2009
:
Eddy transport and mixing in a wind and buoyancy driven jet on the sphere
.
J. Phys. Oceanogr.
,
39
,
1133
1149
.
Cerovečki
,
I.
,
L. D.
Talley
, and
M. R.
Mazloff
,
2011
:
A comparison of Southern Ocean air–sea buoyancy flux from an ocean state estimate with five other products
.
J. Climate.
,
24
,
6283
6306
.
Cessi
,
P.
,
2008
:
An energy-constrained parameterization of eddy buoyancy flux
.
J. Phys. Oceanogr.
,
38
,
1807
1820
.
Cessi
,
P.
,
W. R.
Young
, and
J. A.
Polton
,
2006
:
Control of large-scale heat transport by small-scale mixing
.
J. Phys. Oceanogr.
,
36
,
1877
1895
.
Davis
,
R.
,
1991
:
Observing the general circulation with floats
.
Deep-Sea Res.
,
38A
,
S531
S571
.
de Szoeke
,
R. A.
, and
M. D.
Levine
,
1981
:
The advective flux of heat by mean geostrophic motions in the Southern Ocean
.
Deep-Sea Res.
,
28A
,
1057
1085
.
Döös
,
K.
, and
D.
Webb
,
1994
:
The Deacon cell and other meridional cells of the Southern Ocean
.
J. Phys. Oceanogr.
,
24
,
429
442
.
Eden
,
C.
, and
R.
Greatbatch
,
2008
:
Diapycnal mixing by meso-scale eddies
.
Ocean Modell.
,
23
,
113
120
.
Farneti
,
R.
,
T. L.
Delworth
,
A. J.
Rosati
,
S. M.
Griffies
, and
F.
Zeng
,
2010
:
The role of mesoscale eddies in the rectification of the Southern Ocean response to climate change
.
J. Phys. Oceanogr.
,
40
,
1539
1557
.
Ferrari
,
R.
, and
C.
Wunsch
,
2009
:
Ocean circulation kinetic energy: Reservoirs, sources, and sinks
.
Annu. Rev. Fluid Mech.
,
41
,
253
282
.
Ferrari
,
R.
, and
M.
Nikurashin
,
2010
:
Suppression of eddy diffusivity across jets in the Southern Ocean
.
J. Phys. Oceanogr.
,
40
,
1501
1519
.
Ferreira
,
D.
,
J.
Marshall
, and
P.
Heimbach
,
2005
:
Estimating eddy stresses by fitting dynamics to observations using a residual-mean ocean circulation model and its adjoint
.
J. Phys. Oceanogr.
,
35
,
1891
1910
.
Ganachaud
,
A.
, and
C.
Wunsch
,
2000
:
Improved estimates of global ocean circulation, heat transport and mixing from hydrographic data
.
Nature
,
408
,
453
459
.
Gent
,
P.
, and
J.
McWilliams
,
1990
:
Isopycnal mixing in ocean circulation models
.
J. Phys. Oceanogr.
,
20
,
150
155
.
Gent
,
P.
,
J.
Willebrand
,
T.
McDougal
, and
J.
McWilliams
,
1995
:
Parameterizing eddy-induced tracer transports in ocean circulation models
.
J. Phys. Oceanogr.
,
25
,
463
475
.
Green
,
J. S.
,
1970
:
Transfer properties of the large-scale eddies and the general circulation of the atmosphere
.
Quart. J. Roy. Meteor. Soc.
,
96
,
157
185
.
Hallberg
,
R.
, and
A.
Gnanadesikan
,
2006
:
The role of eddies in determining the structure and response of the wind-driven Southern Hemisphere overturning: Results from the Modeling Eddies in the Southern Ocean (MESO) project
.
J. Phys. Oceanogr.
,
36
,
2232
2252
.
Haney
,
R. L.
,
1971
:
Surface thermal boundary condition for ocean circulation models
.
J. Phys. Oceanogr.
,
1
,
241
248
.
Held
,
I. M.
, and
V. D.
Larichev
,
1996
:
A scaling theory for horizontally homogeneous, baroclinically unstable flow on a beta plane
.
J. Atmos. Sci.
,
53
,
946
953
.
Henning
,
C. C.
, and
G. K.
Vallis
,
2005
:
The effects of mesoscale eddies on the stratification and transport of an ocean with a circumpolar channel
.
J. Phys. Oceanogr.
,
35
,
880
897
.
Hill
,
C.
,
D.
Ferreira
,
J.-M.
Campin
,
J.
Marshall
,
R.
Abernathey
, and
N.
Barrier
,
2011
:
Controlling spurious diapycnal mixing in eddy-resolving height-coordinate ocean models: Insights from virtual deliberate tracer release experiments
.
Ocean Modell.
,
in press
.
Holloway
,
G.
,
1986
:
Estimation of oceanic eddy transports from satellite altimetry
.
Nature
,
323
,
343
344
.
Hughes
,
C. W.
,
1997
:
Comments on “On the obscurantist physics of ‘form drag’ in theorizing about the circumpolar current.”
J. Phys. Oceanogr.
,
27
,
209
210
.
Hutchinson
,
D. K.
,
A. M.
Hogg
, and
J. R.
Blundell
,
2010
:
Southern Ocean response to relative velocity wind stress forcing
.
J. Phys. Oceanogr.
,
40
,
326
340
.
Ito
,
T.
, and
J.
Marshall
,
2008
:
Control of lower-limb overturning circulation in the Southern Ocean by diapycnal mixing and mesoscale eddy transfer
.
J. Phys. Oceanogr.
,
38
,
2832
2845
.
Ivchenko
,
V. O.
,
K. J.
Richards
, and
D. P.
Stevens
,
1996
:
The dynamics of the Antarctic Circumpolar Current
.
J. Phys. Oceanogr.
,
26
,
753
785
.
Johnson
,
G. C.
, and
H. L.
Bryden
,
1989
:
On the size of the Antarctic Circumpolar Current
.
Deep-Sea Res.
,
36
,
39
53
.
Karsten
,
R. H.
, and
J.
Marshall
,
2002a
:
Testing theories of the vertical stratification of the ACC against observations
.
Dyn. Atmos. Oceans
,
36
,
233
246
.
Karsten
,
R. H.
, and
J.
Marshall
,
2002b
:
Constructing the residual circulation of the ACC from observations
.
J. Phys. Oceanogr.
,
32
,
3315
3327
.
Karsten
,
R. H.
,
H.
Jones
, and
J.
Marshall
,
2002
:
The role of eddy transfer in setting the stratification and transport of a circumpolar current
.
J. Phys. Oceanogr.
,
32
,
39
54
.
Keffer
,
T.
, and
G.
Holloway
,
1988
:
Estimating Southern Ocean eddy flux of heat and salt from satellite altimetry
.
Nature
,
332
,
624
626
.
Killworth
,
P. D.
,
1997
:
On the parameterization of eddy transfer: Part I: Theory
.
J. Mar. Res.
,
55
,
1171
1197
.
Kuo
,
A.
,
R. A.
Plumb
, and
J.
Marshall
,
2005
:
Transformed Eulerian-mean theory. Part II: Potential vorticity homogenization and equilibrium of a wind- and buoyancy-driven zonal flow
.
J. Phys. Oceanogr.
,
35
,
175
187
.
Lapeyre
,
G.
, and
I. M.
Held
,
2003
:
Diffusivity, kinetic energy dissipation, and closure theories for the poleward eddy heat flux
.
J. Atmos. Sci.
,
60
,
2907
2917
.
Large
,
W. G.
, and
S. G.
Yeager
,
2009
:
The global climatology of an interannually varying air–sea flux data set
.
Climate Dyn.
,
33
(
2–3
),
341
364
.
Large
,
W. G.
,
J. C.
McWilliams
, and
S. C.
Doney
,
1994
:
Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization
.
Rev. Geophys.
,
32
,
363
403
.
Lumpkin
,
R.
, and
K.
Speer
,
2007
:
Global meridional overturning
.
J. Phys. Oceanogr.
,
37
,
2550
2562
.
Marshall
,
D.
,
1997
:
Subduction of water masses in an eddying ocean
.
J. Mar. Res.
,
55
,
201
222
.
Marshall
,
G.
,
2003
:
Trends in the southern annular mode from observations and reanalyses
.
J. Climate
,
16
,
4134
4144
.
Marshall
,
J.
,
1981
:
On the parameterization of geostrophic eddies in the ocean
.
J. Phys. Oceanogr.
,
11
,
1257
1271
.
Marshall
,
J.
, and
T.
Radko
,
2003
:
Residual mean solutions for the Antarctic Circumpolar Current and its associated overturning circulation
.
J. Phys. Oceanogr.
,
33
,
2341
2354
.
Marshall
,
J.
, and
T.
Radko
,
2006
:
A model of the upper branch of the meridional overturning of the Southern Ocean
.
Prog. Oceanogr.
,
70
,
331
345
.
Marshall
,
J.
, and
K.
Speer
,
2011
:
Closing the meridional overturning circulation through Southern Ocean upwelling
.
Nat. Geosci.
,
in press
.
Marshall
,
J.
,
A.
Adcroft
,
C.
Hill
,
L.
Perelman
, and
C.
Heisey
,
1997a
:
A finite-volume, incompressible Navier Stokes model for studies of the ocean on parallel computers
.
J. Geophys. Res.
,
102
,
5753
5766
.
Marshall
,
J.
,
C.
Hill
,
L.
Perelman
, and
A.
Adcroft
,
1997b
:
Hydrostatic, quasi-hysdrostatic, and non-hydrostatic ocean modeling
.
J. Geophys. Res.
,
102
,
5733
5752
.
Marshall
,
J.
,
E.
Shuckburgh
,
H.
Jones
, and
C.
Hill
,
2006
:
Estimates and implications of surface eddy diffusivity in the Southern Ocean derived from tracer transport
.
J. Phys. Oceanogr.
,
36
,
1806
1821
.
McCarthy
,
M. C.
, and
L. D.
Talley
,
1999
:
Three-dimensional isoneutral potential vorticity structure in the Indian Ocean
.
J. Geophys. Res.
,
104
(
C6
),
13 251
13 267
.
McCartney
,
M. S.
,
1977
:
Subantarctic mode water
.
A Voyage of Discovery, George Deacon 70th Anniversary Volume, M. V. Angel, Ed., Pergamon, 103–119
.
McIntosh
,
P. C.
, and
T. J.
McDougall
,
1996
:
Isopycnal averaging and the residual mean circulations
.
J. Phys. Oceanogr.
,
26
,
1655
1661
.
McWilliams
,
J. C.
,
W. R.
Holland
, and
J. H. S.
Chow
,
1978
:
A description of numerical Antarctic Circumpolar Currents
.
Dyn. Atmos. Oceans
,
2
,
213
291
.
Meredith
,
M. P.
, and
A. M.
Hogg
,
2006
:
Circumpolar response of Southern Ocean eddy activity to a change in the southern annular mode
.
Geophys. Res. Lett.
,
33
,
L16608
,
doi:10.1029/2006GL026499
.
Munk
,
W.
, and
E.
Palmén
,
1951
:
Note on the dynamics of the Antarctic Circumpolar Current
.
Tellus
,
3
,
53
55
.
Naveira Garabato
,
A.
,
R.
Ferrari
, and
K.
Polzin
,
2011
:
Eddy stirring in the Southern Ocean
.
J. Geophys. Res.
,
116
,
C09019
,
doi:10.1029/2010JC006818
.
Olbers
,
D.
,
1998
:
Comments on “On the obscurantist physics of ‘form drag’ in theorizing about the circumpolar current.”
J. Phys. Oceanogr.
,
28
,
1647
1655
.
Olbers
,
D.
,
D.
Borowski
,
C.
Voelker
, and
J.
Wolff
,
2004
:
The dynamical balance, transport and circulation of the Antarctic Circumpolar Current
.
Antarct. Sci.
,
16
,
439
470
.
Plumb
,
R. A.
, and
R.
Ferrari
,
2005
:
Transformed Eulerian-mean theory. Part I: Nonquasigeostrophic theory for eddies on a zonal-mean flow
.
J. Phys. Oceanogr.
,
35
,
165
174
.
Polvani
,
L. M.
,
D. W.
Waugh
,
G. J. P.
Correa
, and
S.-W.
Son
,
2011
:
Stratospheric ozone depletion: The main driver of twentieth-century atmospheric circulation changes in the Southern Hemisphere
.
J. Climate
,
24
,
795
812
.
Prandtl
,
L.
,
1925
:
Bericht untersuchungen zur ausgebildeten turbulenz
.
Z. Angew. Math. Mech.
,
5
,
136
139
.
Prather
,
M. J.
,
1986
:
Numerical advection by conservation of second-order moments
.
J. Geophys. Res.
,
91
(
D6
),
6671
6681
.
Rintoul
,
S.
,
C.
Hughes
, and
D.
Olbers
,
2001
:
The Antarctic Circumpolar Current system
.
Ocean Circulation and Climate, G. Siedler, J. Church, and J. Gould, Eds., Academic Press, 171–302
.
Smith
,
K. S.
, and
J.
Marshall
,
2009
:
Evidence for enhanced eddy mixing at mid-depth in the Southern Ocean
.
J. Phys. Oceanogr.
,
39
,
50
69
.
Speer
,
K.
,
S.
Rintoul
, and
B.
Sloyan
,
2000
:
The diabatic Deacon cell
.
J. Phys. Oceanogr.
,
30
,
3212
3223
.
Stammer
,
D.
,
1998
:
On eddy characteristics, eddy transports, and mean flow properties
.
J. Phys. Oceanogr.
,
28
,
727
739
.
Stone
,
P. H.
,
1972
:
A simplified radiative-dynamical model for the static stability of rotating atmospheres
.
J. Atmos. Sci.
,
29
,
405
417
.
Straub
,
D.
,
1993
:
On the transport and angular momentum balance of channel models of the Antarctic Circumpolar Current
.
J. Phys. Oceanogr.
,
23
,
776
783
.
Talley
,
L. D.
,
2008
:
Freshwater transport estimates and the global overturning circulation: Shallow, deep and throughflow components
.
Prog. Oceanogr.
,
78
,
257
303
.
Taylor
,
G. I.
,
1921
:
Diffusion by continuous movements
.
Proc. London Math. Soc.
,
S2-20
,
196
212
.
Thompson
,
A. F.
,
2008
:
The atmospheric ocean: Eddies and jets in the Antarctic Circumpolar Current
.
Philos. Trans. Roy. Soc. London
,
366A
,
4529
4541
.
Thompson
,
A. F.
, and
W. R.
Young
,
2007
:
Two-layer baroclinic eddy heat fluxes: Zonal flows and energy balance
.
J. Atmos. Sci.
,
64
,
3214
3232
.
Thompson
,
D. W. J.
, and
J. M.
Wallace
,
2000
:
Annular modes in the extratropical circulation. Part II: Trends
.
J. Climate
,
13
,
1018
1036
.
Toggweiler
,
J. R.
,
2009
:
Shifting westerlies
.
Science
,
232
,
1434
1435
.
Toggweiler
,
J. R.
, and
B.
Samuels
,
1998
:
On the ocean’s large-scale circulation near the limit of no vertical mixing
.
J. Phys. Oceanogr.
,
28
,
1832
1853
.
Toggweiler
,
J. R.
, and
J.
Russell
,
2008
:
Ocean circulation in a warming climate
.
Nature
,
451
,
286
288
.
Treguier
,
A. M.
,
I.
Held
, and
V.
Larichev
,
1997
:
Parameterization of quasigeostrophic eddies in primitive equation ocean models
.
J. Phys. Oceanogr.
,
27
,
567
580
.
Treguier
,
A. M.
,
J.
Le Sommer
, and
J. M.
Molines
,
2010
:
Response of the Southern Ocean to the southern annular mode: Interannual variability and multidecadal trend
.
J. Phys. Oceanogr.
,
40
,
1659
1668
.
Viebahn
,
J.
, and
C.
Eden
,
2010
:
Toward the impact of eddies on the response of the Southern Ocean to climate change
.
Ocean Modell.
,
34
,
150
165
.
Visbeck
,
M.
,
J.
Marshall
, and
T.
Haine
,
1997
:
Specification of eddy transfer coefficients in coarse-resolution ocean circulation models
.
J. Phys. Oceanogr.
,
27
,
381
403
.
Watson
,
A. J.
, and
A. C.
Naveira Garabato
,
2006
:
The role of Southern Ocean mixing and upwelling in glacial-interglacial atmospheric CO2 change
.
Tellus
,
58B
,
73
87
.
Wolfe
,
C. L.
, and
P.
Cessi
,
2009
:
Overturning circulation in an eddy-resolving model: The effect of the pole-to-pole temperature gradient
.
J. Phys. Oceanogr.
,
39
,
125
143
.
Wolfe
,
C. L.
, and
P.
Cessi
,
2010
:
What sets the strength of the mid-depth stratification and overturning circulation in eddying ocean models?
J. Phys. Oceanogr.
,
40
,
1520
1538
.
Wolfe
,
C. L.
, and
P.
Cessi
,
2011
:
The adiabatic pole-to-pole overturning circulation
.
J. Phys. Oceanogr.
,
41
,
1795
1810
.

Footnotes

1

This choice of parameters corresponds to a sensitivity of ∂Qeff/∂Ts ~ 15 W m−2 K−1.

2

Ferrari and Nikurashin (2010) recently refined the idea to include the modulation of Le 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.