This study uses a sector configuration of an ocean general circulation model to examine the sensitivity of circumpolar transport and meridional overturning to changes in Southern Ocean wind stress and global diapycnal mixing. At eddy-permitting, and finer, resolution, the sensitivity of circumpolar transport to forcing magnitude is drastically reduced. At sufficiently high resolution, there is little or no sensitivity of circumpolar transport to wind stress, even in the limit of no wind. In contrast, the meridional overturning circulation continues to vary with Southern Ocean wind stress, but with reduced sensitivity in the limit of high wind stress. Both the circumpolar transport and meridional overturning continue to vary with diapycnal diffusivity at all model resolutions. The circumpolar transport becomes less sensitive to changes in diapycnal diffusivity at higher resolution, although sensitivity always remains. In contrast, the overturning circulation is more sensitive to change in diapycnal diffusivity when the resolution is high enough to permit mesoscale eddies.
The Southern Ocean encircles Antarctica and connects the major ocean basins through the agency of the Antarctic Circumpolar Current (ACC) and its associated meridional overturning circulation (MOC). Cold abyssal waters, enriched in carbon and nutrients, upwell in the Southern Ocean amidst a complex interleaving of water masses, giving its circulation a global significance (Rintoul et al. 2001; Meredith et al. 2011). As the crossroads of the oceans, understanding the Southern Ocean circulation, and how that circulation might change, is thus crucial to understanding both the past and future climate of the Earth system.
The simple model due to Gnanadesikan (1999, hereafter G99) heuristically links the global pycnocline depth, and thus the circumpolar transport of the ACC (TACC) through thermal wind balance (Gnanadesikan and Hallberg 2000; Munday et al. 2011), to four processes:
Southern Ocean wind forcing;
the eddy bolus transport, via baroclinic instability, in the Southern Ocean;
deep water formation at Northern high latitudes; and
global diapycnal mixing.
The Southern Hemisphere westerly winds may have been significantly different from the present day mean climate at times in the past (see, for example, Otto-Bliesner et al. 2006). Similarly, estimates of tidal mixing for the Last Glacial Maximum (LGM) suggest that diapycnal mixing was higher (Egbert et al. 2004), particularly in the North Atlantic (Green et al. 2009). However, obtaining robust estimates of global palaeoceanographic circulations, whether at the LGM or otherwise, remains difficult owing to a paucity of data (Wunsch 2003). As a result, numerical and analytical models of varying complexity must be used to assess how such changes might have impacted the Southern Ocean circulation and global climate. Projections of future climate also suggest that changes in both the magnitude and position of the Southern Hemisphere westerlies are plausible (Solomon et al. 2007). The consequences for Southern Ocean circulation and the potential for climate feedbacks have yet to be robustly determined.
In the context of the G99 model, the response of the Southern Ocean circulation to changing forcing has been investigated using general circulation models for wind forcing (Saenko et al. 2002; Fyfe and Saenko 2006; Delworth and Zeng 2008; Allison et al. 2010; Wang et al. 2011), diapycnal diffusivity (Gnanadesikan and Hallberg 2000; Munday et al. 2011), and northern sinking (Fučkar and Vallis 2007). However, these models are usually noneddy resolving, necessitating the use of an eddy parameterization for the entire mesoscale eddy field.
Recent results indicate that the response of an eddy-resolving ocean model to changes in wind forcing differs from that of a noneddy-resolving ocean model with parameterized eddies. In terms of the circumpolar transport, resolving the eddy field leads to a much lower sensitivity to increased wind forcing (Hallberg and Gnanadesikan 2001; Tansley and Marshall 2001; Hallberg and Gnanadesikan 2006). Such lack of sensitivity was first suggested by Straub (1993) and continues to be observed in a growing range of eddy-resolving models (Hogg and Blundell 2006; Meredith and Hogg 2006). This phenomenon has become known as eddy saturation and can be thought of as a marginally critical balance being maintained by the tendency for near-surface Ekman transport to steepen isopycnals and baroclinic eddies to flatten them. Investigations into the eddy saturation behavior of numerical models have recently been extended to primitive equation models using realistically complex geometry (Farneti et al. 2010; Farneti and Delworth 2010). Results indicating a prevailing eddy saturation–type regime for the Southern Ocean and ACC continue to accrue. Several authors (see Hofmann and Morales-Maqueda 2011, for a recent example) have expressed the opinion that eddy saturation must begin at some finite, nonzero wind stress. This has yet to be tested in the limit of zero wind stress.
Consideration of the impact of the southern annular mode (SAM) in an eddy-resolving model also shows a characteristically different response to that of a model with parameterized eddies (Screen et al. 2009). Changes in eddy kinetic energy (EKE) and eddy heat flux take place with a characteristic time scale of 2–3 years, roughly equivalent to that seen in simpler numerical models (Hogg et al. 2008). Furthermore, this time scale is corroborated by that seen in the observational record (Böning et al. 2008; Meredith et al. 2004).
Attention has recently moved to the response of the Southern Ocean MOC to changing wind forcing, which can be broadly decomposed into the wind-driven Eulerian upwelling and the downwelling eddy bolus flux (Johnson and Bryden 1989; Marshall 1997; Marshall and Radko 2003). The net residual MOC (RMOC) is then a subtle balance between these two opposing contributions (Watson and Naveira Garabato 2006). It is the RMOC that is the most relevant circulation when considering, for example, the transport of temperature, nutrients, and other climatically important tracers. Use of the RMOC to describe the upwelling in the Southern Ocean eliminates the spurious “Deacon cell” and its unphysical overestimate of net upwelling (Danabasoglu et al. 1994, also see Döös and Webb 1994).
Hallberg and Gnanadesikan (2006) describe the MOC response to changes in wind forcing in their eddy-permitting hemispherical isopycnal model as “attenuated” with respect to a coarse-resolution version of the same model with parameterized eddies. This reduced response of the RMOC to wind anomalies, when compared to the linear change in the Ekman overturning contribution, was christened eddy compensation by Viebahn and Eden (2010) and, like eddy saturation, has been observed in a growing range of ocean models (Henning and Vallis 2005; Spence et al. 2009; Farneti et al. 2010). In the idealized study of Abernathey et al. (2011) the degree of eddy compensation was found to be crucially dependent upon the form of the surface buoyancy condition, with a fixed flux surface buoyancy forcing resulting in reduced sensitivity to changes in wind stress when compared with a restoring condition. Regardless, in all cases the increase in both upper and lower cells of the RMOC was always much less than the increase in the Eulerian upwelling expected from changes in the surface layer’s Ekman transport.
The response of an ocean general circulation or coupled climate model to changes in the wind stress felt by the ocean component is crucially dependent upon the way the model’s eddy parameterization is formulated. Most noneddy-resolving general circulation models of the ocean and coupled climate models use a derivative of the parameterization due to Gent and McWilliams (1990, hereafter GM). In models with a constant GM coefficient (usually referred to as a thickness diffusivity, but more correctly a quasi-Stokes diffusivity, see Kuhlbrodt et al. 2012), the response to altered wind stress is a strongly correlated change in circumpolar transport and associated changes in stratification and overturning (Fyfe and Saenko 2006). However, a GM-style parameterization formulated such that the GM coefficient can vary with the ocean state, such as that due to Visbeck et al. (1997), breaks this correlation (Kuhlbrodt et al. 2012).
Farneti and Gent (2011) and Gent and Danabasoglu (2011) use different coupled climate models in which the ocean component has a GM-style parameterization where the coefficient can vary in both space and time. This results in an ocean that demonstrates a degree of eddy compensation, such that the response of the overturning to a wind perturbation is reduced with respect to a control run of the same model. Similar results are achieved by Hofmann and Morales-Maqueda (2011) in an ocean-only model using a 3D generalization of Visbeck et al.’s (1997) GM formulation. However, in all three models, the GM parameterization is formulated such that increases in the GM coefficient are tied to steeper isopycnal slopes. Hence, the models may be eddy compensated, but this comes at the expense of a changing circumpolar transport, that is, they are not eddy saturated.
The perturbed coupled climate models of Farneti and Gent (2011) and Gent and Danabasoglu (2011) were run for 100–200 model years, starting from previous experiments. This is short with respect to the spin up time of the ocean, which may extend to multimillennial time-scales (Wunsch and Heimbach 2008; Allison et al. 2011; Jones et al. 2011). The ocean-only model of Hofmann and Morales-Maqueda (2011) was run for substantially longer (a 3000-yr control run with 1000-yr perturbation experiments), which should place it closer to its final equilibrium state. It is difficult to assess how the circumpolar transport in the perturbed models might vary at statistical steady state. However, based on the increase in circumpolar transport seen by Hofmann and Morales-Maqueda (2011) [a 35-Sv (1 Sv ≡ 106 m3 s−1) increase for a doubling of wind stress] and the continued upward trajectory of the transport in Fig. 9 of Farneti and Gent (2011) for their experiments with increased wind stress, it seems unlikely that eddy saturation would eventually occur. The extent to which these phenomena occur concurrently in an eddying ocean is therefore crucial in guiding future developments of eddy parameterizations.
Eddy saturation and eddy compensation are not the same phenomenon (Meredith et al. 2012), and the connection between them remains the subject of ongoing research. However, taken together they present an exciting possibility regarding climate change. The close link between stratification and transport, highlighted by the model of G99, indicates that eddy saturation is not just an argument regarding circumpolar transport; an eddy-saturated circumpolar current would have very close to invariant isopycnal slope/stratification. A totally eddy-compensated ocean would also show little change in up-/downwelling along those isopycnals. Hence, modification of global climate by increased upwelling of carbon-rich water in the Southern Ocean, as hypothesized by Toggweiler et al. (2006), might well have only limited applicability to the real climate system.
The wind stress perturbations applied in the experiments discussed above have typically been modest, and so the full range of parameter space has not been sampled. Furthermore, to capture all the time scales inherent to ocean adjustment, the full meridional extent of the ocean must be simulated (Allison et al. 2011), and the model integrated to equilibrium to enable a robust understanding of the results (Treguier et al. 2010). Current ocean-only and coupled climate models are run with a wide range of diapycnal diffusivities, and little is really known regarding the sensitivity of an eddy-resolving model to this parameter. While the degree of sensitivity to the globally integrated diapycnal diffusivity is crucially dependent upon the presence or absence of wind (Munday et al. 2011), assessment of changes in this sensitivity in eddy-resolving models has yet to be carried out over the same swath of parameter space. Here we investigate both of these problems in an eddy-permitting ocean-only configuration of the Massachusetts Institute of Technology (MIT) general circulation model (MITgcm) that has been integrated to as close to equilibrium as we can afford.
We begin in section 2 by describing our computational domain and some of the numerical choices and compromises made. The key results from a wide range of numerical experiments are introduced in section 3. A more detailed discussion of the experiments begins with a summary of the basic model state at a range of grid spacings in section 4. Sections 5 and 6 discuss two separate sets of experiments in which the surface wind stress and diapycnal diffusivity have been independently varied. We close with a summary of our conclusions and a discussion of their significance in section 7.
2. Model numerics and domain
To investigate the effects of eddy saturation and eddy compensation at equilibrium, we employ the MITgcm (Marshall et al. 1997a,b) in an idealized modeling framework. This allows us to conduct many experiments under different combinations of forcing and parameter choices. Geometrically, the model domain takes the form of the narrow sector shown in Fig. 1 in spherical polar coordinates, with the implied full variation of the Coriolis parameter. Latitudinally this domain extends from 60°S to 60°N, but is only 20° wide in longitude. At the southern end of the sector, a re-entrant channel of 20° latitudinal width allows for the formation of a circumpolar current. The use of such model geometry has a considerable provenance within the ocean-modeling literature as a computationally efficient means to investigate, for example, the meridional overturning circulation (Bryan 1986; von der Heydt and Dijkstra 2008; Wolfe and Cessi 2009), the parameter sensitivity of general circulation models (Bryan 1987), the processes setting ocean stratification (Henning and Vallis 2004, 2005; Wolfe and Cessi 2010), and the impact of ocean circulation on ocean carbon storage and atmospheric carbon dioxide concentration (Ito and Follows 2003, 2005). Key model parameters, which are discussed in more detail below, are given in Table 1.
The model domain’s meridional extent is dictated by the desire to produce an interhemispheric overturning of the type currently thought to take place in the North Atlantic. This prevents a vanishingly small residual overturning in the deep ocean caused by a solid boundary at the equator. It also avoids the use of artificial deepening of the stratification on the northern boundary of the circumpolar channel using, for example, northern boundary restoring zones (Hallberg and Gnanadesikan 2006; Viebahn and Eden 2010; Abernathey et al. 2011; Morrison et al. 2011) or increased diapycnal diffusivity (Hogg 2010). Extending the domain beyond the equator includes important adjustment time scales, which might otherwise be excluded, and allows for deep-water formation at the northern boundary and diapycnal mixing to play an important role in the volume budget of the global pycnocline in a G99-style budget (Allison et al. 2011). It comes at an enhanced computational cost because of the increase in the size of the domain; this disadvantage is more than offset by the ability of the stratification and RMOC to evolve together dynamically.
Achieving equilibrium is important in understanding the role of eddies in ocean models (Treguier et al. 2010), but remains elusive in costly global model simulations. Our decision to use a longitudinally narrow sector is dictated by the desire to integrate experiments with a fine enough grid spacing to permit/resolve eddies to thermodynamic and biogeochemical equilibrium. By restricting the domain in this way, to roughly a factor of 10–20 times smaller than the actual longitudinal extent of the ocean, we are able to sufficiently speed up our model runs to attain this goal. Early tests with a narrower sector indicate that 20° is sufficiently wide to prevent undue influence on the eddy field of the domain width.
The meridional extent of the circumpolar channel is several times that of Drake Passage, although it remains narrower than the width of the Southern Ocean south of Australia and New Zealand. Initial low-resolution testing demonstrated that a channel of this width was required to obtain a circumpolar transport of the same order as that observed for the ACC (not shown). Similarly wide circumpolar channels are often used in other sector model applications (see, for example, Henning and Vallis 2005; Ito and Follows 2005).
The bathymetry of the model is idealized to a flat-bottom over the majority of the domain. Throughout the main 20° wide sector, the ocean has a depth of 5000 m. However, a 1° (or 1 grid box, whatever is the greater) region of half depth is appended to the domain within the confines of the circumpolar channel (see Fig. 1). This ridge blocks all contours of f/H, where f is the Coriolis frequency, and H is the depth, and ensures that the zonal momentum input by the wind is balanced by bottom form stresses (Munk and Palmén 1951; Johnson and Bryden 1989).
In the vertical, the model domain has 42 unevenly spaced levels. The surface layer is the thinnest, at 10 m, while the bottom 10 levels have a common thickness of 250 m (the maximum used). All model resolutions use a constant grid spacing, in degrees, in the zonal direction. However, to ensure accurate representation of the eddy field, the grid spacing in the meridional direction is scaled by the cosine of latitude. This means that, regardless of latitude, the grid boxes are approximately square.
All model resolutions use the Gent–McWilliams (GM) parameterization of mesoscale eddies (GM; Gent et al. 1995) via the skew-flux formulation of Griffies (1998), which also includes the rotation of “vertical” diffusivity into the diapycnal direction (Redi 1982). At coarse grid spacing, Δ = 2° or coarser, GM is used as a parameterization of the entire eddy field with the coefficient (κGM) set to a large value of 1000 m2 s−1, typical of that used in general circulation models.1 At a grid spacing of Δ = ½° and finer, GM is used as an adiabatic subgridscale closure for the turbulent dissipation of potential temperature and salinity, as per the argument of Roberts and Marshall (1998). The coefficient value is chosen such that the contribution of GM to the meridional heat transport is small with respect to the other terms. See Table 2 for the values of κGM used at different grid-spacings. No other forms of buoyancy dissipation are otherwise needed. Near the surface the parameterization is tapered as per Gerdes et al. (1991) with a maximum isopycnal slope of 2 × 10−3 applied. This combination was chosen to prevent the shutdown of high-latitude convection via GM “taking over” in the near surface layers at coarse resolution.
We use the vector invariant form of the momentum equations with the wet point method of Jamart and Ozer (1986) applied to the Coriolis terms. Potential temperature and salinity are advected using the one-step monotonicity-preserving seventh-order advection scheme due to Daru and Tenaud (2004). This reduces numerical diapycnal mixing to a low level (Hill et al. 2012). Many tests were carried out at eddy-permitting grid spacings to choose the most appropriate tracer advection scheme. These tests had to be run for several hundred model years to give a reliable indication of the level of numerical diapycnal diffusivity. Side and bottom boundary conditions are no slip. Convective adjustment is carried out by increasing the vertical diffusivity to 100 m2 s−1 whenever static instability is present. Viscous dissipation closures for the momentum equations are resolution dependent, with harmonic friction being used for grid spacings of 2° and coarser, and biharmonic friction at grid spacings of ½° and ⅙°. Regardless of grid spacing, viscosity is specified via a grid-scale Reynolds number, which ensures that as resolution and time step are covaried, the viscosity also changes appropriately. As with the choice of tracer advection scheme, appropriate specification of this parameter is essential in minimizing numerical diapycnal mixing (Ilicak et al. 2012). As with the choice of tracer advection scheme, many tests had to be performed to choose the most appropriate viscosity. The model uses the nonlinear equation of state of Jackett and McDougall (1995) with temporally constant reference pressure. Model parameters that are independent of grid spacing are given in Table 1, with resolution-dependent model parameters in Table 2.
In keeping with the idealized geometry, we choose similarly idealized surface forcing and boundary conditions. Wind forcing is applied as an idealized jet with a peak value of 0.2 N m−2 at 45°S. The shape of this jet is given by
where θ is the latitude, and τ0 is the peak wind stress. The wind stress field used in the basic state experiments is shown graphically in Fig. 2a. There is no wind significantly north of the reentrant channel to emphasize the role of the Southern Ocean winds in ventilating the deep ocean and driving Southern Ocean circulation. Because of the narrowness of the sector, gyres generated by wind to the north would be severely reduced in transport with respect to, for example, the observed North Atlantic subtropical gyre. Furthermore, the use of wind stress to the north could potentially introduce a basin contribution to the circumpolar transport, as identified by Nadeau and Straub (2009, 2012). For these reasons, we choose to exclude such gyres from the forcing entirely. This wind forcing is very similar to that used by Viebahn and Eden (2010). The key difference is the presence of some wind stress curl just to the north of the circumpolar channel, which helps connect the overturning within the channel to that of the basin. The off-center position of the peak wind forcing, with respect to the center of the channel, places the model ACC toward the northern edge of the channel. There is considerable precedent for the use of such abbreviated wind stress, including Viebahn and Eden (2010), Morrison et al. (2011), Morrison and Hogg (2013), and Nikurashin and Vallis (2011, 2012). Furthermore, in experiments with coupled models in which the wind stress is externally altered, such as those due to Farneti et al. (2010), it is typically the wind stress only in the circumpolar belt that is substantially modified.
For potential temperature (henceforth temperature), we choose a simple restoring condition with an idealized profile similar to that used in many multiple equilibria studies (see, for example, von der Heydt and Dijkstra 2008). This condition broadly mimics the observed surface temperature distribution of the Atlantic, with cold water at the northern/southern boundary and warm water at the equator. The functional form of the surface restoring temperature is given by
where TS is the temperature at the southern boundary, ΔT is the southern boundary-to-equator temperature difference, and TN is the temperature at the northern boundary. For the experiments reported here, we choose TS = 0°C, ΔT = 30°C, and TN = 5°C, which is broadly similar to the observed temperature distribution of the Atlantic (although with an exaggerated surface temperature at the equator). This temperature restoring condition is illustrated in Fig. 2b. This basic state forcing is comparable to the warm pole experiment of Wolfe and Cessi (2009).
As with temperature, we take a simple restoring approach to the surface salinity of our model ocean. The condition is chosen to give increased salinity in the tropics, broadly consistent with the observed surface distribution of salinity. By placing the highest salinity at the equator, rather than displaced to the north and south as in the North Atlantic, we ensure that salt is exported from the Tropical regions. The profile of restoring salinity is given by
where SS is the salinity at the southern boundary, ΔS is the southern boundary-to-equator salinity difference, and SN is the salinity at the northern boundary. For the experiments reported here, we take SS = 34, ΔS = 3, and SN = 34, so that the salinity at the northern and southern boundary is the same. The meridional profile of salinity restoring condition is shown in Fig. 2c. More correctly we should prescribe the surface freshwater flux, but by retaining a restoring condition, we allow ourselves the luxury of being able to determine the surface density, and thus the stratification and overturning regime, of the sector model a priori (see below for details.). We also reduce the possibility of a salt advection feedback operating on the meridional overturning, leading to multiple equilibria (Stommel 1961; Rahmstorf and Ganopolski 1999; Johnson et al. 2007).
The combination of the above surface restoring conditions for temperature and salinity gives rise to the surface-referenced potential density (henceforth density) restoring profile illustrated in Fig. 2d. These profiles have been chosen such that the restoring density at the northern boundary lies between the restoring density at the southern boundary and that at the northern edge of the reentrant channel (as shown by the dotted line in Fig. 2d). The relative densities of the water masses formed at these three latitudes then results in southern boundary water forming an abyssal water mass that serves as a model analog of Antarctic Bottom Water (AABW). Arrayed vertically above this abyssal water mass, northern boundary water forms an analog of North Atlantic Deep Water (NADW), while water from the northern edge of the re-entrant channel forms Antarctic Intermediate Water (AAIW) at the shallowest depth. The overturning associated with this stratification looks like the North Atlantic; there is deep water formation at both northern and southern boundaries and two cells, much like those described by Toggweiler et al. (2006) (see section 4c).
By perturbing the surface density restoring profile, typically by setting SN to values higher or lower than SS, it is possible to produce two other broad categories of stratification. A sufficient increase in SN will create the densest water in the domain at the northern boundary. This shuts down the bottom (AABW) cell of the RMOC and produces a warm, salty abyssal water mass almost solely made up of model analog NADW. Water from the southern end of the domain becomes surface trapped, forming a cold, fresh surface layer within the circumpolar channel. In contrast, if SN is selected such that the water at the northern boundary is less dense than any of the water within the channel, the model analog NADW becomes surface trapped; the entire domain becomes dominated by model analog AABW and AAIW.
Regardless of the surface-restoring conditions, it is possible to obtain qualitatively similar RMOCs/stratifications to these alternative states by perturbing the model parameters. For example shutting off the wind or significantly increasing κυ will both produce the second type of RMOC/stratification described above (see sections 5b and 6b, respectively). The RMOCs/stratifications described here are much like those of Vallis (2000), and Henning and Vallis (2004, 2005), and indicate that the gross parameters of the ocean circulation that select for the particular oceanic RMOC/stratification remain the same regardless of whether the simulated ocean is eddying in nature.
In the sections that follow we discuss several suites of model experiments that use grid spacings of 2°, ½°, and ⅙°. The coarsest of these grid spacings is noneddy resolving, while at ½° and ⅙° closed vortices and filaments of vorticity/fronts in temperature and salinity are generated, respectively (see section 4a for details). The ⅙° grid does not resolve scales finer than the first Rossby deformation radius. However, in this case the grid is fine enough to resolve the eddies themselves, if not the scale at which they form at, and fronts and filaments of vorticity have started to appear in the model results (see section 4). At each grid spacing, two suites of experiments are discussed. In the first, the peak wind stress (τ0) is varied, around the basic state value of 0.2 N m−2, from 0 to 1 N m−2. In the second suite, the diapycnal diffusivity is varied over the range 1 − 30 × 10−5 m2 s−1, a range typical of that used in noneddy-resolving models for a variety of applications. Regardless of model resolution, the basic state experiments use a diapycnal diffusivity of 3 × 10−5 m2 s−1. To shorten the integration time required, each suite was initialized from the final stratification of a previous, lower resolution set of experiments interpolated to the new grid spacing. The 2° were initialized from a set of very coarse 4° experiments and the ½° experiments were then initialized from the result of the 2° experiments. After 1000 years, the ½° results were then interpolated to ⅙°, and these experiments begun.2 Where time-average results are discussed, the 2° experiments have been averaged over 1000 years, the ½° over 100 years, and the ⅙° over 10 years.
3. Key results
The key results of our numerical experiments are summarized in Fig. 3, where the relationship between the time-mean “circumpolar” transport (the zonal transport through the re-entrant channel) and the strength of the wind forcing (Fig. 3a) and diapycnal diffusivity (Fig. 3b) are shown. Different averaging periods are used for each grid spacing; 1000 years for 2°, 100 years for ½°, and 10 years for ⅙°. The bars represent two standard deviations of the instantaneous monthly transport about the mean. They indicate the instantaneous variability of the circumpolar current, rather than the standard error in the mean, which is extremely small due to the large number of sample values in the averaging period.
Examination of Fig. 3a demonstrates that the noneddy-resolving model (2°, blue line) behaves like other global climate models employing a constant GM coefficient, that is, the circumpolar transport changes strongly with the wind stress (Fyfe and Saenko 2006). Even with no wind at all (τ0 = 0 N m−2) a significant TACC of ~50 Sv occurs. This transport occurs for the reasons elucidated by Munday et al. (2011), that is, that the pycnocline to the north of the ACC is deepened by diapycnal mixing, even in the absence of wind. This then leads to a considerable circumpolar transport via thermal wind shear. The increase in TACC with wind forcing continues across the extreme range considered here, which reaches a peak wind stress of 1.0 N m−2, compared to the basic state value of 0.2 N m−2. The increase in transport does not remain linear with wind stress, although it is close to this limit across many of the experiments. The reader should note that no error bars are shown on the Δ = 2° line of Fig. 3a as the variability is so low that they would be smaller than the plotted symbol in most cases.
When the grid spacing is refined to ½° (red line), and again to ⅙° (green line), the model behaves like the high-resolution numerical models discussed in section 1. In other words, TACC “saturates” at some finite value of wind stress and ceases to increase with further increases in wind stress. Indeed, for the first time our ⅙° experiments demonstrate that such saturation may take place with no wind at all, since the increase in variability effectively makes the green line on Fig. 3a indistinguishable from flat. The extreme range of wind forcing considered in the experiments presented here also demonstrates that an increase in wind stress cannot overwhelm the eddy processes responsible for the eddy saturation phenomenon and that it can be expected to be present even outside of the observed wind stress range. Further refinement of the model grid spacing seems unlikely to lead to a change in this conclusion.
The circumpolar transport of the zero wind stress case at both ½° and ⅙° is substantial; for both of these grid spacings it is ~2 × the transport at 2°. This is due to the way a resolved, or at least permitted, eddy field can respond to changes in wind forcing. At 2°, the GM parameterization is as capable of flattening the isopycnals at zero wind as it is at 5× the wind, since the coefficient is constant. At ½° and ⅙°, the eddy field weakens as the wind stress decreases, and so the isopycnals are not flattened quite as much. Hence, the thermal wind across the channel, which is fixed because of the restoring condition, and the transport of the circumpolar current, is much higher than in the 2° zero wind case.
The increase in variability with increasing wind for the ½° and ⅙° experiments in Fig. 3a shows that as the wind stress increases, the extra energy from the increased surface wind work is transferred into the variability of the system, rather than the mean state. This makes it difficult to assess the saturation state of a given circumpolar current by examining instantaneous snapshots. For example, at the strongest wind forcing (equivalent to a peak wind stress of 1.0 N m−2) at a grid spacing of ⅙° the circumpolar transport may instantaneously reach ~200 Sv, or be as low as ~0 Sv, even though the long-term average transport is ~100 Sv. This occurs due to the movement of eddies across the virtual “mooring” of the line of latitude used to calculate the transport. When taking into account the extreme variability present in the experiments with higher wind stress, it is easy to appreciate that the resulting trend line may have an instantaneous gradient spanning a wide range of magnitudes and be of either sign. As a result, monitoring of the ACC to determine its saturation state may require many years of records, to generate a robust long-term average and to further determine whether this average is changing. Although the extreme variability for the strongest wind forcing is well outside of the estimated variability of the ACC (see, for example, Meredith et al. 2011, and references therein), we note that Nadeau and Straub (2009) also see a large increase in variability of circumpolar transport in a two-layer quasigeostrophic model. In the case of low bottom drag, Nadeau and Straub (2012) attribute this to the presence of a basin mode in their solution and large-scale barotropic Rossby waves breaking at the western boundary. In the quasigeostrophic results of Nadeau and Straub (2009, 2012), the circumpolar transport is not eddy saturated from the zero wind limit, unlike in our experiments.
Altering the diapycnal diffusivity used by the model causes TACC to vary as shown in Fig. 3b. The range of diapycnal diffusivity used in these experiments is similar to that used in a range of general circulation models for a variety of purposes, such as ocean-only investigations of the ocean’s general circulation and the response of a coupled ocean–atmosphere climate model to increased concentrations of greenhouse gases. Over this range of diffusivity, the model gives a significant variation in mean circumpolar transport for all three model resolutions. In contrast to the wind experiments, the 2° and ½° sets of experiments exhibit approximately the same sensitivity to changes in diapycnal diffusivity. There is increased variability in the ½° experiments, relative to the 2° experiments, owing to the generation of vortices by undamped instability and an accompanying increase in internal variability.
As with the wind experiments, the model exhibits a clear separation between the behavior of the mean TACC at ⅙° and that at 2° and ½°. There is also an increase in the variability of the instantaneous circumpolar transport, although, it is not sufficient to argue that the ⅙° line is effectively flat. However, the gradient of the line is clearly reduced with respect to the other grid spacings. This begs the question: if the grid spacing was sufficiently refined, would the model demonstrate complete insensitivity of the mean circumpolar transport in response to changes in diapycnal mixing?
In summary, at the highest resolution used, a grid spacing of ⅙°, the model shows different behavior in response to perturbations in wind stress and diapycnal diffusivity. Under changing wind forcing, it demonstrates complete eddy saturation from zero wind stress upward. In contrast, when the diapycnal diffusivity is varied, there is still significant variation in the mean circumpolar transport. We emphasize that eddy saturation is a steady state/time-mean argument; at ½° and ⅙° grid spacing, the model always exhibits intense variability in the instantaneous TACC. Furthermore, the extent to which the circumpolar currents are eddy compensated is also dependent upon both the type and magnitude of the applied forcing perturbation. See sections 5 and 6 for further details.
4. Overview of the basic state
a. Mean flow and stratification
The transition from noneddy-resolving (2°), through ½°, to eddy-permitting (⅙°) grid spacings alters both the instantaneous and long-term time-average flow. At 2° grid spacing, the flow is viscous and dissipative. The eddy kinetic energy (EKE, shown in Fig. 4a) is not due to the propagation of mesoscale eddies, but rather to large-scale, low-frequency variations of the geostrophic flow. The result is an EKE field at least four orders of magnitude smaller than that of the higher model resolutions used. This variability does increase with both wind forcing and diapycnal diffusivity, but always remains about 104 times smaller than the ½° EKE. The circumpolar current is very broad and consists of a single “jet” extending over 5°–10° of latitude. Inspection of individual snapshots of the barotropic streamfunction shows that they are all much like the 1000-yr average shown in Fig. 4a, with no closed contours of streamfunction forming owing to the absence of strong vortices. Similarly, the instantaneous temperature field always resembles the surface temperature shown in Fig. 4d, once a statistical steady state has been reached.
At both ½° and ⅙°, the EKE is truly due to the presence of a vigorous mesoscale eddy field. The EKE increases significantly between these two grid spacings, both in terms of magnitude and spatial extent, as the ⅙° grid can resolve more of the mesoscale eddy field. At ⅙° the viscous dissipation of momentum is now sufficiently low for the steepness of the isopycnals near the northern boundary to lead to undamped instability processes and the generation of a vigorous eddy field, as shown by the surface temperature field in Figs. 4e and 4f. Examination of the instantaneous temperature and salinity fields also reveals that the very low strength western boundary currents, formed via a thermal wind response to the meridional variation of stratification and resulting inflow/outflow into the western boundary region at low/high latitudes, also become unstable. While here the focus is on Southern Ocean/southern boundary processes, these regions could potentially lead to important eddy fluxes of both buoyancy and biogeochemical tracers elsewhere in the domain. At ⅙°, Fig. 4f also demonstrates the increase in variability at low latitudes, as Rossby waves are continuously excited at the eastern boundary and propagate westward.
The contours of barotropic streamfunction in Figs. 4b and 4c also reveal the qualitative change in form that the barotropic flow undergoes as the grid spacing is refined. Even at ½°, the circumpolar flow is concentrated into a much narrower jet across the northern edge of the channel. At ⅙°, the jet is narrower still. The small region of negative streamfunction to the north of the channel is due to the curl of the wind stress creating a very small gyre. At the eddying resolutions, the mean jet also shows the presence of a standing wave pattern, which becomes stronger for the finer grid. By ⅙° the standing wave is strong enough to create closed contours of time-mean streamfunction. The instantaneous barotropic streamfunction at ½° and ⅙° shows considerable deviation from the time-mean picture of Figs. 4b and 4c, since the circumpolar current itself can migrate and the eddies propagate. In addition to general instability throughout the channel, the standing wave undergoes instability processes and barotropic eddies, visible as closed contours of barotropic streamfunction, are able to propagate southward along the edge of the step within the channel (not shown). The instantaneous surface temperature in the ½° and ⅙° simulations, shown in Figs. 4e and 4f, respectively, show many of these features. In particular, Fig. 4e indicates that large-scale closed vortices are generated at ½°. At ⅙°, Fig. 4f indicates the generation of fronts in the temperature field, which can be seen as filaments in the vorticity field (not shown).
It is notable in Fig. 3 that there is a systematic reduction in TACC between the ⅙° experiments and the coarser grid spacings. This is due to a reduction in the thermal wind shear, which ultimately determines much of the circumpolar transport due to the very weak bottom flow. This is dominated by changes in global pycnocline depth, although a slight warming of the abyss also contributes. The two isopycnals drawn in Fig. 5 for all three basic states indicate that while the light density layers are at approximately the same depth for each grid spacing, the depth of the heavier water masses can vary considerably between grid spacings. In particular, we see that the effect of the strong surface restoring is to pin the density outcrops to roughly the same latitude, but that in the ⅙° experiment the 1036.5 kg m−3 isopycnal is shallower at the northern edge of the channel. The reduced slope of this isopycnal across the channel and thus smaller meridional density gradients, gives rise to a weaker thermal wind shear and lower circumpolar transport for the ⅙° basic state, and the other ⅙° experiments (see Fig. 3a). The layering of the water masses that form the vertical stratification is as described at the end of section 2, that is, analogs of AABW and AAIW form at the southern boundary and the northern edge of the re-entrant channel, with an analog of NADW intruding between them from the northern boundary. The use of a restoring condition for salinity does limit the strength of any comparison with the observed North Atlantic. Notably, it prevents the NADW from having high salinity, and so this water mass does not appear as a salinity maximum in the vertical stratification.
In the absence of any continental barriers, the vertically and zonally integrated (indicated by chevrons) zonal momentum budget of the ocean can be written as (Johnson and Bryden 1989; Stevens and Ivchenko 1997)
where τ is the surface wind stress, 〈τb〉 is the (total) bottom stress, F(ub) is the bottom drag, pb is the bottom pressure, and is the bottom slope.3 The last term on the right-hand side is the form stress caused by the variation of pressure across bathymetry.
In our model, the presence of the shelf in the reentrant channel blocks all of the f/H contours. This prevents the model from developing a strong bottom flow and results in TACC being determined almost purely by the stratification (Wolff et al. 1991; Marshall 1995). As a result, 〈F(ub)〉 in Eq. (4) is small and the prevailing momentum balance has the zonal momentum input from the wind balanced by form stress across the shelf. This is a fundamentally different balance to that in Viebahn and Eden (2010) and Abernathey et al. (2011), whose flat-bottom allows for the development of strong bottom flow and an absence of bottom form stress.
b. Adjustment time scales
Treguier et al. (2010) make the point that understanding an eddying ocean’s response to SAM wind events requires that these models be at thermodynamic equilibrium. Reaching such an equilibrium is a major motivation for the use of an idealized geometry such as ours. To assess any remaining drift in the model experiments, we have examined a number of metrics, including domain average temperature/salinity and the circumpolar transport. Typically, there is no question as to the steady state nature of the 2° simulations, as they can be trivially run for as long as required, in this case 20 000 years. However, determination of the steady state nature, or otherwise, of the ½° and ⅙° experiments is, at least partly, subjective. For this reason, Figs. 6a and 6b present the domain average potential temperature as a function of time for the ½° and ⅙° basic states, respectively, as well as the extreme wind stress experiments.
The ½° experiments have all been run for 4000 years. As Fig. 6a demonstrates, the domain average temperature shows little meaningful drift by this point. For the more weakly forced experiments, here characterized by the zero wind stress case, the required spin up time is significantly extended. Indeed, after 1000 years, the zero wind stress experiment is clearly still adjusting. This long centennial-to-millennial adjustment time scale is predicted by conceptual models (Allison et al. 2011; Jones et al. 2011) and is a consequence of the time it takes for a small change in the net volume transport out of the Southern Ocean to alter the properties of the large volume of water to the north of the circumpolar latitudes. In the case of the ⅙° experiments, the difference between the initial and final model state is smaller than for both 2° and ½°, that is, a smaller modification is required. Hence, there is a noticeable decrease in the spin-up time as resolution increases.
The spinup of the ⅙° experiments is characterized in Fig. 6b. Again, the more strongly forced experiments show considerably faster adjustment, at least during the initial phase. After 400 years, the point at which the results here are shown, the domain-average temperature is still varying. However, we have continued the three experiments in Fig. 6b, and the experiments with the highest and lowest diapycnal diffusivities (30 × 10−5 m2 s−1 and 1 × 10−5 m2 s−1, respectively), for a further 400 years, for a total of 800 years of model integration. By this point the domain average temperature is typically drifting by approximately ±0.01°C century−1, or less, with the experiment with five times the basic state wind drifting at about twice this rate.
Of more relevance to our main conclusions (see section 3) is the degree to which the mean circumpolar transport is still varying. The instantaneous variation of this diagnostic, as a series of monthly values, is shown for the century/decade averaging period of the ½° and ⅙° experiments in Figs. 7a and 7b, respectively. This demonstrates the clear separation between the different experiments at ½°, but the overlapping nature of TACC with different wind stresses at ⅙°. Table 3 indicates that over the second 400 years of the extended model run, the mean TACC does change in all 5 of the experiments. However, it is only the zero wind stress experiment that does not have its mean TACC after 400 years within one standard deviation of the mean TACC after 800 years. The extended run of the ⅙° basic state and 4 extreme experiments indicate that our conclusions regarding the straightness of the ⅙° line in Fig. 3a and the curvature present in Fig. 3b are not altered by the low degree of transience still present in the stratification.
c. Residual meridional overturning circulation
In an eddying ocean the circulation that actively transports tracers is not simply the Eulerian-mean of the velocity field. Rather, it includes a “bolus” component, which can be derived from the thickness-weighted mean flow (Gent et al. 1995; Lee et al. 1997) and is equivalent to the residual velocity in a transformed Eulerian-mean formulation (McIntosh and McDougall 1996; Marshall 1997). Care must be taken as to the method by which the residual circulation is calculated, with quantitatively different circulations resulting from different methods of averaging (Nurser and Lee 2004a,b).
The RMOC is here diagnosed from the model results using the thickness-weighted averaging method of Abernathey et al. (2011), that is, the following integral is calculated at every time step
where h = ∂z/∂ρ is the layer thickness in potential density referenced to the 30th model level (at ~2000 m depth, henceforth “potential density”). The integral is calculated using 241 discrete layers that are 0.025 kg m−3 thick. It is this streamfunction that describes the circulation of tracers most accurately, eliminating the “Deacon cell” in the Southern Ocean (Danabasoglu et al. 1994) and closely resembling the structure of the circulation in an isopycnal model (see, for example, Hallberg and Gnanadesikan 2006). Using this method of calculation combines the Eulerian-mean circulation caused by surface Ekman transport and eddy-induced bolus flux into a single diagnostic.
Broadly speaking, the residual overturning for each of the different model resolutions is similar. In all cases, a model analog of the upper NADW cell overlies a lower AABW cell. As seen in Figs. 8a–c, the upper cell does not significantly change topology or strength with finer grid spacing. The potential density at which water flows into the reentrant channel does lighten slightly at ⅙°, but is otherwise relatively invariant, given the large differences in EKE between the experiments. The upper cell leads to upwelling within the confines of the re-entrant channel (shown by the dashed line), with subduction actively occurring at the northern edge. This gives rise to a model analog of the Antarctic Intermediate Water (AAIW) cell, which is very much surface trapped. At ⅙° this corresponds to a local peak in the convective index (not shown), demonstrating shallow, local convection driving this third cell to form mode waters. Notably, the bottom cell strengthens with increasing model resolution. We attribute this to the increase in strength of the southward eddy flux of heat across the geostrophically blocked band of latitudes. This gives rise to a large increase in surface heat flux, due to the restoring condition on temperature, and deep convection, which takes place within the last few grid boxes next to the southern boundary, and thus drives a stronger circulation. In general, the bottom cell remains weak, due to the absence of significant deep water formation to the south of “Drake Passage,” in this model setup.
The thick gray contour in Figs. 8a–c is the time- and zonal-average surface potential density. This contour shows little variation between resolutions because of the salinity and temperature restoring conditions applied at the surface. The presence of streamlines of the RMOC above the time and zonal mean potential density contour in Figs. 8a–c is an indication that, in an instantaneous sense, important diabatic transformations take place at densities lower than the time-zonal-mean. In areas where eddy processes are strong, such as in the reentrant channel, these transformations are stronger and/or more frequent, resulting in more streamlines above the gray contour. Near the northern boundary the advection of potential density anomalies anticlockwise around the basin can also lead to this effect.
5. Sensitivity to wind forcing
As the numerical grid spacing is decreased in the narrow sector model, the circumpolar transport eventually becomes invariant to changes in wind stress (see section 3). However, other aspects of the ocean circulation, such as the residual overturning circulation or the statistical properties of the eddy field, may remain sensitive to changes in wind forcing. Here we examine such aspects of the model’s circulation for the experiments in which the wind stress is varied, that is, those used in the production of Fig. 3a.
a. Eddy kinetic energy
In Fig. 4 it is clear that as the model grid spacing is refined, EKE increases, owing to improved resolution of mesoscale instability processes and the eddies they generate. This is symptomatic of the more vigorous mesoscale eddy field and is common in most general circulation models of the ocean as the eddy-resolving threshold is passed (see, for example, Roberts et al. 2004; Farneti et al. 2010). Typically this increase in EKE is seen as an improvement in realism, with respect to observations of EKE in the real ocean by, for example, satellite altimetry. As Fig. 9a demonstrates, increasing the wind, at a fixed grid spacing, leads to an increase in EKE. The large separation between the EKE of equivalent ½° and ⅙° experiments remains, although the EKE of the ½° experiment with five times the basic state wind (peak wind stress of 1.0 N m−2 compared to 0.2 N m−2) has comparable EKE to the ⅙° basic state. The EKE of 2° experiments always remains at least three orders of magnitude lower than the ½° experiments, and so is not shown. For the ⅙° experiments, the growth of the EKE with the wind forcing is slightly faster than linear.
Figure 9a demonstrates that once closed vortices, and internal variability, begin to occur at ½°, there is a strong correlation between surface wind forcing and the domain average EKE. When combined with Fig. 3a, it shows that the increasing rate of wind work at the surface acts to increase the eddy kinetic energy, rather than the mean kinetic energy. This is unsurprising, given that the dominant mechanical energy budget, for a zonal channel in the adiabatic limit of weak time-mean bottom flow, is given by (Cessi et al. 2006; Cessi 2008; Abernathey et al. 2011)
where is the time-average surface flow, is the bottom EKE, and G indicates a function of the bottom EKE, vertical viscosity, etc, consistent with the bottom no-slip boundary condition. The growth of EKE with wind forcing is then a consequence of the fact that the surface wind-work against the time-mean current must be balanced by mechanical dissipation acting on the bottom EKE. This is consistent with the quasi-linear nature of the dependence of EKE on the wind stress. Diagnostics confirm that this balance holds to leading order (not shown); indeed, the average EKE in the bottom layer is little different in magnitude to the values shown in Fig. 9a.
The time-mean EKE does not allow us to distinguish between larger eddies and stronger eddies, although inspection of sea surface height and temperature fields indicates that the eddies do indeed become both larger and stronger with stronger wind stress. The movement of these eddies then leads to the much larger variation in circumpolar current observed in Fig. 3a, and the stronger vertical momentum/southward heat transport required to close the form drag balance of the circumpolar current (Olbers 1998).
b. Residual meridional overturning circulation
Viebahn and Eden (2010) find that the sensitivity of the overturning to changes in wind forcing in their model is strongly dependent upon the size of the grid spacing, that is, the presence of a resolved eddy field. Typically, at higher resolution (finer grid spacing), the sensitivity is reduced. However, even at a 5-km grid spacing, their model does not produce total eddy compensation, that is, the RMOC does still change in response to an altered wind field. The eddy-resolving channel model of Abernathey et al. (2011) demonstrates that sensitivity of the residual overturning circulation to changes in wind forcing crucially depends on the formulation of the surface buoyancy forcing condition. In particular, the sensitivity is enhanced for relaxation boundary conditions, as used by ourselves and Viebahn and Eden (2010), with respect to fixed flux boundary conditions. This is due to changes in the implied surface buoyancy flux under a restoring boundary condition, which the residual circulation must match (Walin 1982). As a result, such modification of the fluxes has a significant impact upon the RMOC. The sensitivity of the overturning to surface buoyancy forcing was recently investigated by Morrison et al. (2011). As with Viebahn and Eden (2010), sensitivity to changes in the magnitude of the surface flux of buoyancy was reduced at finer grid spacing. However, these three studies have all considered domains that were truncated to the north, requiring that the northern boundary stratification be prescribed in some way. In contrast, here the stratification to the north is allowed to freely vary.
The RMOC for the lowest (no wind) and highest (peak wind stress of 1.0 N m−2) perturbation experiments are shown in Fig. 10 for all three grid spacings. The RMOCs have been calculated in the same way as those for the basic state, as described in section 4c. There is a large qualitative change in the type of circulation that results in the absence of wind (see Figs. 10a–c); the quasi-adiabatic inflow of the NADW cell no longer penetrates into the reentrant channel, which results in a circulation dominated by the combined AABW and AAIW cells. This is due to the lack of Ekman transport in the surface model layer. Removal of this transport removes the suction on the lower layers that results in the NADW cell penetrating into the channel. As such, the eddy bolus transport, which does not go to zero owing to finite isopycnal slope and eddy generation, defines the resulting residual circulation. This regime change is independent of model resolution, as one might expect, although quantitative differences between the resolutions still occur.
When the surface wind stress is increased, the residual overturning tends to strengthen, as shown in Figs. 10d–f. The RMOC is still quasi-adiabatic, outside of the surface layers, but the inflow into the channel and subsequent upwelling has increased. This general increase occurs regardless of model resolution. However, visual inspection indicates that the RMOC of the ½° and ⅙° experiments (Figs. 10e and 10f, respectively) has increased to a lesser degree than the 2° experiment. As with the basic state RMOC of section 4c, the bottom/AABW cells noticeably increase in strength for these extreme wind forcing simulations. This is due to the increase in convection at the southern boundary, rather than any direct wind-driving effects. However, the integrated zonal momentum balance of the current remains of the form in Eq. (4).
In Fig. 9b the change in the RMOC as surface wind stress is altered is quantified by taking the maximum streamfunction value north of the equator as a measure of the upper (NADW) cell and the minimum between 30°S and the equator as a measure of the lower (AABW) cell. In both calculations the surface layers, which may become noisy particularly at 2° grid spacing, are avoided. Qualitatively speaking, these two measures show broadly the same picture as other measures, for example the streamfunction at an appropriately chosen potential density and latitude.
Figure 9b indicates that regardless of model resolution, the upper cell retains finite sensitivity to wind forcing. However, the ½° and ⅙° lines have a shallower gradient then the 2° line regardless of the wind stress. This is particularly noticeable for the ⅙° experiments. This points toward the possibility of total eddy compensation in the residual overturning, just as total eddy saturation occurred for the circumpolar volume transport in section 3. Viebahn and Eden (2010) and Morrison and Hogg (2013) both see a similar loss of sensitivity of the RMOC to wind stress changes in their model experiments, although not quite to the same degree. This could be for a combination of reasons, including the cropping of their model domains at the equator (both), the use of a linear equation of state (both), the use of fixed flux surface buoyancy forcing (Morrison and Hogg 2013), the fundamentally different momentum balance (Viebahn and Eden 2010), and the more moderate wind forcing perturbations they apply (both). The wind forcing perturbations considered here extend well beyond the range the cited papers consider.
In contrast to the upper cell, the lower cell (dashed lines in Fig. 9b) begins to show increased sensitivity to the wind stress at high wind stress multiples for the ⅙° experiments. Sensitivity of the bottom cell remains low in the 2° and ½° experiments. We attribute this to the same mechanism that gives rise to a slight increase in the bottom cell in the basic state as the model grid is refined, that is, the southward eddy heat flux increases, warms the water at the southern boundary, and results in a feedback on the restoring condition to give a stronger heat flux. The sign of the change in the overturning of the bottom cell is in the opposite sense to that of Abernathey et al. (2011, see their Fig. 5), who see an increasingly positive streamfunction in their bottom cell as wind increases under restoring boundary conditions. This could be due to the large changes in experimental design, specifically the use by Abernathey et al. of a meridionally truncated domain with full-depth restoring zone on the northern boundary. Anecdotally, we have found that such restoring at depth is considerably more effective at preventing changes in the abyssal stratification, when compared to a surface restoring condition. As a result, our abyssal stratification can change to a larger degree between both model resolutions and individual experiments at the same grid spacing. Such changes could lead a considerably different residual flow, potentially allowing for the increase in the bottom cell that we see here.
c. Pycnocline depth and stratification
where D is the pycnocline depth, σ2 is the potential density referenced to the 30th model level (the same reference level used for calculation of the RMOC), and z = −H is the depth of the lowest model level. This is essentially a center-of-mass calculation for the vertical coordinate. After area-weighted averaging between 30°S and 30°N, this provides a quantitative measure of how deep the pycnocline is. Averaging over a limited area of the domain in this way avoids the northern and southern boundaries, where a deep mixed layer may develop and bias the resulting diagnostic heavily toward these regions. The result of this calculation for the wind experiments at all three grid spacings is shown in Fig. 9c.
Figure 9c summarizes the contents of Fig. 5, in as much as its shows that the pycnocline depth in the ⅙° basic state, characterized by the σ2 = 36.5 kg m−3 isopycnal in Fig. 5, is shallower than both the 2° and ½° pycnocline. It also confirms that the pycnocline depth in the 2° and ½° basic state calculations are effectively the same. Primarily, the figure reinforces the close link between pycnocline depth and circumpolar transport highlighted in section 1, that is, that through thermal wind they are strongly coupled and, thus, the eddy saturation hypothesis is as much about pycnocline depth/stratification as it is about transport. It indicates that while the pycnocline depth deepens considerably at 2° under increased wind forcing, the extent of deepening at both ½° and ⅙° is much reduced. For example, at 2° the difference between pycnocline depth at zero wind and 5 multiples of the basic state wind is ~1800 m. However, at ⅙° the difference between these most extreme pycnocline depths is only ~150 m.
6. Sensitivity to diapycnal diffusivity
As the numerical grid spacing is decreased in the narrow sector model, the circumpolar transport becomes less sensitive to changes in diapycnal diffusivity (see section 3). In contrast to the wind experiments, the model does not reach a point of total insensitivity to further changes in diapycnal diffusivity at a grid spacing of ⅙°. Here we examine other aspects of the ocean circulation already discussed for the wind experiments in section 5, for experiments in which the diapycnal diffusivity has been both increased and decreased.
a. Eddy kinetic energy
The variation of EKE with diapycnal diffusivity coefficient is shown in Fig. 11a (as with Fig. 9a, the 2° experiments are omitted for clarity due to the near-vanishing EKE). As with increased wind forcing, an increase in the coefficient gives rise to an increase in EKE, although not to the same extent. It also remains true that the more refined grid of the ⅙° experiments is more able to represent the mesoscale eddy field and the instability that generates it. As a result, the EKE for the ⅙° experiments is always higher than for the equivalent ½° experiment. Similarly, the western boundary currents themselves are increasingly unstable, as is the flow close to the northern boundary.
With reference to the conceptual model of G99, an increased diapycnal diffusivity acts to deepen the pycnocline. Making the reasonable assumption that the pycnocline will outcrop within the confines of the channel, this will lead to a steepening of the isopycnals across the channel. Hence, the flow will become more susceptible to instability and higher EKE will result. At the grid spacings considered here, the increase in resolution is not sufficient to completely offset the increased diapycnal diffusivity. As a result, the pycnocline deepens (see section 6c), as the increased diffusivity acts as a source of potential energy.
In terms of the approximate mechanical energy budget given by Eq. (6), both the surface wind work, due to increased surface current speed, and dissipation of bottom EKE increase with κυ. A second term should also be added to the right-hand side of Eq. (6) to account for the rate at which diapycnal mixing converts potential energy to kinetic energy.
b. Residual meridional overturning circulation
The RMOCs that result from altering the diapycnal diffusivity coefficient are shown for the extreme values of κυ = 1 × 10−5 m2 s−1 and κυ = 30 × 10−5 m2 s−1 at all three model resolutions in Fig. 12. As with the altered wind experiments, we find a qualitative change in the regime of the overturning from an “Atlantic” type to a “Pacific” type, that is, penetration of the NADW cell into the channel region ceases and the AAIW and AABW cells expand and combine together. However, this regime change is in the opposite sense to that brought about by wind forcing; as the “forcing” increases (κυ increases) the upper cell becomes surface trapped and the lower cell dominates the circulation. In the case of the wind experiment, as “forcing” decreases (wind stress decreases) a qualitatively similar regime change occurred. While the change in regime in superficially similar to that with zero wind, the rate of overturning increases as κυ increases, just as it would with increased wind stress. Furthermore, due to the way that diapycnal diffusivity affects the flow, there is no way to increase the rate of overturning while maintaining an adiabatic regime, as occurs with increased wind forcing.
Reduction of the diapycnal diffusivity leads to an increasingly adiabatic flow, see Figs. 12a–c for the κυ = 1 × 10−5 m2 s−1 simulations, since the prescribed flow across isopycnals is now reduced. This remains true regardless of the model grid spacing. At the coarsest grid spacing used, 2°, numerical problems begin to occur, in which abyssal water masses can form that are colder than the coldest surface water. This leads to a section of the domain that is isolated from the surface, since convection cannot penetrate into this cold, dense water mass. The net effect is that seen in Fig. 12a, that the bottom cell is substantially reduced in size and magnitude. Care was taken to avoid this numerical issue in all other model runs, although it was found to be unavoidable at this particular combination of low κυ and coarse model grid.
Quantification of the overturning in Fig. 11b, using the maximum/minimum of the RMOC streamfunction for the upper/lower (NADW/AABW) cell as described in section 5b, reveals that the model resolutions tend to diverge from each other at high diapycnal diffusivities. The strength of the upper cell tends to increase with κυ, even as it becomes restricted in spatial extent at high diapycnal diffusivity. As with the wind experiments, the upper cell is less sensitive to changes in diapycnal diffusivity at ⅙° grid spacing. This indicates that eddy processes continue to buffer the NADW cell against changes even when the cell no longer extends into the reentrant channel.
Figure 11b also reveals that the bottom cell is indeed strengthening. In fact, it shows that the ⅙° experiments are actually more sensitive to changes in κυ than the coarser grid spacings. Essentially, the increase in diapycnal mixing is able to supply more potential energy to both cells, however, it does so to the lower cell preferentially, such that this cell becomes dominant. The change in overturning that results from changes in potential energy sources/sinks is dependent upon the spatial distribution of those sources and sinks (Oliver and Edwards 2008). As a result, it is nontrivial to determine the changes in overturning that might result from changes in diapycnal diffusivity within the Southern Ocean, as opposed to the global change in diffusivity used here.
In contrast to the case for wind changes, the AABW cell of the RMOC is not compensated with regards to diapycnal diffusivity changes, as shown in Figs. 12 and 11b. In fact, Fig. 11b indicates exactly the opposite to eddy compensation behavior, since an increase in model resolution and EKE leads to a steeper slope, and thus increased sensitivity to changes in diapycnal diffusivity. This is because the Eulerian and bolus components of the RMOC are no longer in competition. Instead, the eddy bolus overturning reinforces the Eulerian overturning. The increase in the bolus overturning due to changes in resolution/model parameter/forcing can then act to increase the RMOC, rather than oppose the accompanying change in the Eulerian overturning. Hence, the RMOC becomes stronger in an eddy-resolving calculation once the overturning has switched to this diabatic regime.
c. Pycnocline depth and stratification
As revealed by Fig. 11c, the range of pycnocline depth achieved at all three grid spacings is broadly comparable; the dramatic reduction in range seen between 2° and ½° for the wind experiments (Fig. 9c) is not evident. Across all of the values of κυ considered, the pycnocline depth in the 2° experiments varies by over 300 m, as opposed to ~2700 m for the wind experiments. For both ½° and ⅙°, this decreases slightly; the change in stratification and transport are far less dramatic than the change in the RMOC. Thus the changes in stratification and transport are both consistent with arguments based on thermal wind, that is, that the changes to the transport must be reflected by, or reflect, changes in the stratification.
Recent eddy-resolving numerical simulations of the ocean indicate that the circumpolar transport of the Antarctic Circumpolar Current (ACC), or equivalently the global pycnocline depth, may be remarkably insensitive to changes in wind forcing. Similarly, the residual overturning of the Southern Ocean may also show a certain degree of insensitivity to wind changes, or at least be less sensitive than its wind-driven Eulerian upwelling component. These two phenomena have been christened eddy saturation and eddy compensation, respectively. The two phenomena may occur separately, and the interrelationship between the two remains a topic of research.
Using the MITgcm in an eddy-permitting configuration, we have investigated eddy saturation and eddy compensation under changes in wind forcing and the sensitivity of circumpolar transport/residual circulation to changes in diapycnal diffusivity. Our numerical simulations are, by necessity, idealized. However, there are several key differences to previously published work, such as that of Viebahn and Eden (2010), Abernathey et al. (2011), and Morrison et al. (2011). First, the model domain extends across the equator, such that the low-latitude stratification is not specified a priori. Second, the applied forcing/parameter perturbations span a wide swath of parameter space, such that the robustness and asymptotic behavior of the phenomena can be assessed. Third, the equation of state of the ocean is nonlinear, such that we have made the first steps toward a full treatment of the impact of salinity on the ACC and Southern Ocean circulation.
At sufficiently high resolution, in this case an eddy-permitting grid spacing of ⅙°, the model’s time-mean circumpolar transport shows (almost) complete insensitivity to wind forcing. Crucially, we have shown that this can occur in the limit of zero wind forcing, implying that an ocean with a vigorous mesoscale eddy field may show no variation of its long-term time-mean ACC transport regardless of the presence of wind. We emphasize that eddy saturation is very much an argument about the time-mean state of the Southern Ocean. Even though the time-mean circumpolar transport is almost invariant with changes to wind forcing, we find significant changes to the variability as wind is altered.
As one would expect from thermal wind balance, for experiments in which the time-mean zonal transport shows little variation with wind, the global pycnocline depth is also roughly constant. Instantaneously, the mesoscale eddy field may cause individual isopycnals to heave up-and-down by hundreds of meters. However, away from the circumpolar channel, that is, at low latitudes, the position of the pycnocline does not vary significantly.
In general, higher-resolution models show stronger eddy compensation, with the degree of compensation that takes place increasing with wind stress. Finite sensitivity would always be expected to remain at “low” wind, since the eddy bolus component of the overturning will not go to zero at zero wind, even though the wind driven upwelling would. However, this implies that there may be an upper limit to this finite sensitivity.
Contrary to the results of Abernathey et al. (2011, see their Fig. 5, red circles) using a restoring surface boundary condition, we find that the Antarctic Bottom Water Cell (AABW)/downwelling cell becomes stronger with increased wind, that is, the streamfunction at a particular point in latitude and density becomes more negative. We attribute this to the increasingly strong eddy heat flux near the southern boundary, which leads to increased southern boundary cooling. While this is not directly wind-driven, since the Ekman transport is northward at every latitude and the details may well be model dependent, it may be a secondary effect of the changing wind forcing.
In the case of the set of experiments in which the diapycnal diffusivity is modified from the basic state value, a rather different set of conclusions result. We find that for the three grid spacings considered here, changes to diapycnal diffusivity still result in significant changes in the circumpolar transport and global pycnocline depth. However, there is a weakening of the sensitivity of the transport at finer grid spacing. This suggests that a further refinement in grid spacing might result in the achievement of insensitivity to diapycnal diffusivity of the zonal transport and stratification. The finest grid spacing used here is an eddy-permitting ⅙°, and so we suggest that passing the poorly-defined eddy-resolving threshold might be necessary. Based upon a deformation radius of ~8–10 km in the Southern Ocean, this would require a model grid spacing of no coarser than ~4–5 km, roughly equivalent to .
In terms of the residual overturning, the diapycnal diffusivity experiments demonstrate that an increase in model resolution leads to an increase in sensitivity to parameter changes. This is just the opposite of that expected under eddy compensation for wind forcing changes. In this case, the transition to a diabatic overturning regime creates a situation, within the channel region, in which the Eulerian and bolus components are both acting in the same direction. As a result, there is no chance of compensation between them. Rather, a refinement in the model grid, and thus an increase in EKE and the vigour of the mesoscale eddy field, leads to an increase in the bolus component, and thus in the total residual overturning. From the perspective of which parameter value to use, this indicates that it might be even more crucial for an eddy-resolving model to use the “right” diapycnal diffusivity than one with parameterized eddies.
With reference to the contrast between the wind and diapycnal diffusivity experiments, there may be a global versus local forcing dichotomy. The wind forcing is very localized to within the channel, whereas diapycnal diffusivity also functions throughout the basin. It remains an open question as to whether a localized region of high/low diffusivity within the confines of the channel might have the same effect as the same magnitude of spatially localized variation at the northern boundary. This is particularly relevant given the spatially and temporally varying nature of diapycnal mixing that recent observations have revealed (Polzin et al. 1997; Kunze et al. 2006; Damerell et al. 2012). Such observations indicate that the Southern Ocean is a region of widespread turbulent mixing (Naveira Garabato et al. 2004). Given the trend we see when moving from ½° to ⅙°, it might be possible for a sufficiently high resolution numerical model to become insensitive to changes in diapycnal diffusivity if the increased level of diapycnal diffusivity was only spatially localized to a few specific latitudes, for example.
There are many limitations to the idealized approach that we have taken here. For example, the inclusion of salinity and use of a nonlinear equation of state is certainly a step forwards. However, we have adopted a restoring condition for salinity, which incorrectly allows the local salinity to impact the virtual freshwater flux and eliminates the salinity maximum in the model’s NADW. In our case, this was done to aid the a priori determination of the stratification and overturning. However, a rich range of multiple equilibria can exist in an ocean model with a pure flux condition for salinity. Furthermore, the geometry has been simplified to a single basin. Whether this has effectively condensed the whole global ocean into a 20° sector, or whether it simply represents a narrow Atlantic Ocean is debatable. As with the use of flux forcing for salinity, the use of multiple basins increases the range of available multiple equilibria and would give further insight into how eddy saturation and eddy compensation relate to the complexity of the real world. Ideally, the model resolution would be higher, as with Morrison and Hogg (2013), and the model would be global. However, this would likely then prevent our goal of quasi-equilibria being reach.
The use of a carefully designed eddy parameterization can allow a coarse-resolution climate/ocean general circulation model to reproduce aspects of the eddy-resolving result of eddy compensation, for example, Hofmann and Morales-Maqueda (2011), Farneti and Gent (2011), and Gent and Danabasoglu (2011). However, even these superior schemes calculate the eddy transfer coefficient as a function of stratification. Thus, while eddy compensation can be achieved, it comes at the expense of changing stratification and circumpolar transport, that is, eddy saturation is intrinsically neglected. The use of functions of higher powers of the stratification to specify the coefficients may mitigate the loss of eddy saturation somewhat (Jones et al. 2011). However, such schemes become increasingly arbitrary. Connecting the coefficients used by parameterizations directly to the prognostic model fields, as opposed to arbitrarily varying them with, for example, wind stress, is clearly a desirable end. As suggested by Farneti and Gent (2011), it may be profitable to pursue schemes such as that due to Eden and Greatbatch (2008) and Marshall and Adcroft (2010), which carries EKE as a model variable and allow production terms to be tied directly to wind stress.
The conceptual model of the ocean’s RMOC and stratification due to Nikurashin and Vallis (2012) produces parametric scalings inline with those of the algebraic model of G99. In the limit of weak diapycnal mixing/strong wind, and neglecting the eddy contribution to net Southern Ocean upwelling, the scaling for the overturning predicts that it varies , while the depth scale of the stratification is predict to vary as . Similar scalings for the strong diapycnal mixing/weak wind limit suggests that the overturning scales as , with the depth-scale of the stratification going as . Nikurashin and Vallis (2012) support these scalings with the results of two different numerical models at coarse grid spacing.
Our numerical results indicate that an eddying model is not bound by the scaling rules of Nikurashin and Vallis (2012). The wind stress experiments indicate that beyond the eddy-resolving threshold, the stratification/circumpolar transport would more appropriately scale as ~τ0 and that the overturning would also scale with a low power of the wind stress. This is superficially similar to the weak wind limit result of Nikurashin and Vallis (2012). The diapycnal diffusivity experiments demonstrate that these aspects of the circulation continue to vary with other parameter choices. This implies that a single scaling rule is unlikely to describe the behavior of both coarse and eddy-resolving models. The results of our wind forcing experiments are actually a much better validation of the simple idea of Straub (1993); that beyond a threshold value of TACC, baroclinicity of the ocean prevents any further sensitivity to wind forcing.
In deriving their model, Nikurashin and Vallis (2012) invoke a constant eddy diffusivity, to parameterize the transport effects of the eddy field. Clearly, such an assumption does not hold once eddies have begun to appear in the flow; Abernathey et al. (2011) use mixing length arguments to relate EKE and eddy diffusivity such that they both increase with surface wind stress. Assuming that defining such a diffusivity is realistic, then the variation of EKE with both wind stress and diapycnal diffusivity (Figs. 9a and 11a, respectively) demonstrates that the resulting diffusivity would vary with the forcing/parameter choice and model grid spacing. Hence, we have made no attempt to fit a particular power law to our results, since we expect such a fit to change following subsequent refinement of the model grid. Until numerical general circulation models have reached a resolution high enough to demonstrate that their solutions are converged, we anticipate that such variations of power law will continue. Even in the presence of a robust theory for what sets such eddy diffusivities, numerical convergence would still be required for mutual validation to be achieved.
The climatic implications of the combination of eddy saturation and eddy compensation are important. A totally eddy-saturated circumpolar current would have close to invariant mean zonal transport and stratification, that is, isopycnal slope and depth. If such a current was also totally eddy compensated, the up-/downwelling along these isopycnals would also vary only weakly. Hence, changing wind could do little to increase the rate of upwelling, and subsequent outgassing of abyssal carbon reservoirs through upwelling in the Southern Ocean. The combination of these phenomena therefore places doubt upon the extent to which changes in Southern Ocean wind forcing could alter global climate. However, climatic insensitivity to Southern Ocean wind forcing is contingent upon the earth currently being in a “high” wind regime, such that the Ekman transport does not fall off rapidly with decreasing wind, while the eddy bolus transport stays finitely high. The increased sensitivity of our eddying model results to changes in diapycnal diffusivity suggest that changing such mixing might be a more reliable way in which to alter climate. However, good estimates of diapycnal diffusivity throughout the global ocean, and more information on how an eddy-resolving model responds to localized changes in diapycnal diffusivity, are required to quantitatively assess whether such changes could, for example, help explain the glacial/interglacial change in atmospheric carbon dioxide concentration.
Stuart Daines aided with selection of the most appropriate tracer advection scheme. DRM acknowledges Ryan Abernathey for his MITgcm layers package, which made calculation of the residual overturning quick and simple. Discussions with Andy Hogg, Adele Morrison, Mehmet Ilicak, Raf Ferrari, and Ric Williams are gratefully acknowledged. This work is supported by NERC. HLJ is supported by a Royal Society University Research Fellowship. DPM acknowledges additional support from the Oxford Martin School. This work made use of the facilities of HECToR, the UK’s national high-performance computing service, which is provided by UoE HPCx Ltd at the University of Edinburgh, Cray Inc and NAG Ltd, and funded by the Office of Science and Technology through EPSRC’s High End Computing Programme.
For reasons of numerical stability it was found to be easier to initialize the ⅙° diapycnal diffusivity experiments from the 4° experiments used to initialize the 2° experiments. In some cases, this leads to a noticeable lag between the ⅙° basic state and the 12 experiments that make up the rest of the ⅙° diapycnal diffusivity suite.
The model configuration used here has a no-slip bottom boundary condition, rather than an applied linear or quadratic drag. Hence, the bottom drag term is a function of the vertical viscosity, bottom velocity, etc, that is conveniently summarized in this way.