The impact of Southern Ocean residual upwelling on atmospheric CO2 on centennial and millennial timescales

The Southern Ocean plays a pivotal role in climate change by exchanging heat and carbon, and provides the primary window for the global deep ocean to communicate with the atmosphere. There has been a widespread focus on explaining atmospheric CO2 changes in terms of changes in wind forcing in the Southern Ocean. Here, we develop a dynamically-motivated metric, the residual upwelling, that measures the primary effect of Southern Ocean dynamics on atmospheric CO2 on centennial to millennial timescales by determining the communication with the deep ocean. The metric encapsulates the combined, net effect of winds and air–sea buoyancy forcing on both the upper and lower overturning cells, which have been invoked as explaining atmospheric CO2 changes for the present day and glacial-interglacial changes. The skill of the metric is assessed by employing suites of idealized ocean model experiments, including parameterized and explicitly simulated eddies, with online biogeochemistry and integrated for 10,000 years to equilibrium. Increased residual upwelling drives elevated atmospheric CO2 at a rate of typically 1–1.5 parts per million/106 m3 s−1 by enhancing the communication between the atmosphere and deep ocean. This metric can be used to interpret the long-term effect of Southern Ocean dynamics on the natural carbon cycle and atmospheric CO2, alongside other metrics, such as involving the proportion of preformed nutrients and the extent of sea ice cover.


Introduction
The Southern Ocean plays a fundamental role in the climate system, ventilating much of the global ocean by forming subtropical mode, intermediate and bottom waters, as well as uniquely returning deep waters to the surface (Marshall and Speer 2012). By providing the window for the global deep ocean to communicate with the atmosphere, these different physical transport pathways lead to a delicate carbon balance over the Southern Ocean : atmospheric CO 2 can be drawn down by subduction of upper ocean waters, as well as increased by outgassing of CO 2 from natural carbon-rich waters.
Changes in physical forcing over the Southern Ocean have been invoked as explaining changes in atmospheric CO 2 for the past (Sigman and Boyle 2001;Toggweiler et al. 2006;Menviel et al. 2008;Tschumi et al. 2008Tschumi et al. , 2011d'Orgeville et al. 2010;Sigman et al. 2010) and the present day Morrison et al. 2015;Zickfeld et al. 2007). For example, stronger winds have been invoked in the present day as driving a weakening in the sink of anthropogenic carbon, inducing enhanced surface Abstract The Southern Ocean plays a pivotal role in climate change by exchanging heat and carbon, and provides the primary window for the global deep ocean to communicate with the atmosphere. There has been a widespread focus on explaining atmospheric CO 2 changes in terms of changes in wind forcing in the Southern Ocean. Here, we develop a dynamically-motivated metric, the residual upwelling, that measures the primary effect of Southern Ocean dynamics on atmospheric CO 2 on centennial to millennial timescales by determining the communication with the deep ocean. The metric encapsulates the combined, net effect of winds and air-sea buoyancy forcing on both the upper and lower overturning cells, which have been invoked as explaining atmospheric CO 2 changes for the present day and glacial-interglacial changes. The skill of the metric is assessed by employing suites of idealized ocean model experiments, including parameterized and explicitly simulated eddies, with online biogeochemistry and integrated for 10,000 years to equilibrium. Increased residual upwelling drives elevated atmospheric CO 2 at mixing (Le Quéré et al. 2007), increased upwelling (Lenton andMatear 2007) and a stronger Eulerian-mean overturning (Lovenduski et al. 2007(Lovenduski et al. , 2008, which increases the exposure of natural carbon-rich deep waters to the surface, opposing anthropogenic carbon uptake. Changes in the zonal symmetry of Southern Hemisphere winds have recently been implicated in the reinvigoration of the Southern Ocean carbon sink (Landschützer et al. 2015). Stronger winds have also been proposed as strengthening upwelling and initiating the deglacial rise in atmospheric CO 2 (Anderson et al. 2009). However, the paleoclimatic inferences of implied upwelling have often been interpreted in terms of wind-induced changes in the Eulerian-mean circulation, rather than accounting for the partial compensation from the ocean eddy circulation.
There is a debate over the expected response of the Southern Ocean circulation to altered wind forcing. From diagnostics of high-resolution numerical models, strengthening winds may lead to little overall increase in strength of the large-scale Southern Ocean horizontal circulation, usually inferred from zonal Antarctic Circumpolar Current (ACC) transports, where the wind energy instead increases eddy kinetic energy at the mesoscale (e.g. Gnanadesikan 2001, 2006;Meredith and Hogg 2006); this process is referred to as eddy saturation. An increase in wind forcing is also expected to increase the Eulerian-mean overturning, which may be compensated by an increase in the opposing eddy-driven overturning (Böning et al. 2008;Viebahn and Eden 2010;Abernathey et al. 2011;Wolfe and Cessi 2011); this process is referred to as eddy compensation. Eddy saturation and eddy compensation need not always be tightly coupled. For example, some high-resolution model integrations reveal the eddy saturated limit is approached, but not the eddy compensation limit, so that increasing Southern Ocean winds may enhance the rate of overturning without large changes in ACC transports (Meredith et al. 2012;Morrison and Hogg 2013;Munday et al. 2013). Furthermore, the response of Southern Ocean largescale circulation and its density structure to changes in physical forcing is likely to occur on centennial to millennial timescales. Consequently, any changes in circulation are unlikely to be detected in the relatively short available time series of observations or high-resolution model integrations (Allison et al. 2011;Jones et al. 2011;Shakespeare and Hogg 2012). The lack of a unified view of how physical forcing affects atmospheric CO 2 on longer timescales is hampering our understanding of how to place a context for past and ongoing climate variability.
Motivated by the importance of the Southern Ocean for climate change in the present day and the past, here we propose a new metric-the strength of Southern Ocean residual upwelling-to represent the dominant effect of Southern Ocean dynamics on atmospheric CO 2 on centennial to millennial timescales. Here our focus is on changes in the natural carbon cycle, so that ocean carbon inventory changes equal opposing atmospheric CO 2 changes. The residual upwelling is defined by the upwelling over the Southern Ocean, including the effects of both the upper and lower overturning cells (Fig. 1), accounting for the combined effects of the Eulerian-mean and eddy circulations. The upper cell consists of Upper Circumpolar Deep Water drawn to the surface, which is then transported northward across the ACC and then subducted as Antarctic Intermediate and Subantarctic Mode Water (Talley 2013). This transformation to mode waters requires a buoyancy input from heating or freshwater fluxes to lighten the upwelled waters. The lower cell involves upwelling confined south of the ACC, transferring Lower Circumpolar Deep Water to the surface, which then becomes denser and contributes to the formation of bottom waters adjacent to Antarctica (Talley 2013). This transformation of surface to bottom waters requires a buoyancy loss from cooling and brine rejection, which may be sensitive to the strength of the polar easterly winds (Stewart and Thompson 2012). Our metric of residual upwelling is designed to measure the vertical transfer of dense waters to the upper ocean, and thus is a measure of Fig. 1 A schematic view of the residual upwelling connecting the deep ocean to the atmosphere. The residual upwelling involves both an upper cell, called the residual circulation, including the combined effect of a northward Ekman upper transport and a southward eddy upper transport, and a lower cell, including sinking and northward spreading of bottom waters. There is a CO 2 outgassing to the atmosphere (vertical arrows) from the upwelling of carbon rich-deep waters and a CO 2 influx into the ocean through the subduction of upper ocean waters and formation of bottom waters. The thick full arrows represent overturning streamlines and the thin lines the potential density surfaces. The subtropical upper waters (light red) are warm and contain low carbon and nutrients, while the colder deeper and bottom waters (white and blue) are nutrient and carbon rich the communication of the deep ocean and the atmosphere over the Southern Ocean. This definition differs from alternative views focusing on parts of the overturning circulation, either the wind-induced Eulerian mean upwelling or the residual circulation, neither of which fully account for the combined effects of the upper and lower overturning cells.
To understand our metric of residual upwelling, consider the effect of stronger winds, expected to strengthen the northward upper ocean transport and associated overturning carried by the Eulerian-mean circulation. Stronger winds act to steepen the tilt of northward-deepening density surfaces across the Southern Ocean (Fig. 1). However, in response, ocean eddies act to flatten density surfaces, inhibiting this steepening, and provide an opposing southward upper ocean transport (Marshall 1997;Marshall and Radko 2006;Viebahn and Eden 2010). This partial compensation leads to the residual circulation, defined by the sum of the Eulerian-mean and eddy circulations, only weakly increasing with strengthening winds (Munday et al. 2013). This residual circulation provides the overturning transport of heat and carbon across the ACC, connecting the Southern Ocean with the northern ocean basins. However, this residual circulation does not represent the entire upwelling, as there are upwelling cells in the Southern Ocean that do not extend across the ACC. Accordingly, our metric, the residual upwelling, defining the total upward transport over the Southern Ocean, is made up of the sum of the residual circulation across the Southern Ocean and upwelling cells contained south of the ACC. The residual upwelling then provides a primary measure of the extent that deep waters are upwelled in the open ocean and allowed to communicate with the atmosphere before returning to the ocean interior. The residual upwelling leads to outgassing of CO 2 where carbon-rich deep waters are exposed to the atmosphere and, conversely, ocean uptake of atmospheric CO 2 in watermass formation regions. Our metric captures the balance of these opposing air-sea flux contributions as the system approaches a long-term equilibrium.
In this study, the dynamical control of the open ocean on atmospheric CO 2 is investigated using two different models with contrasting domains (Sect. 2). Both models are integrated for 10,000 years to reach thermodynamic and biogeochemical equilibrium for a range of different surface wind stresses and associated buoyancy forcings. The relationship between residual upwelling and atmospheric CO 2 is investigated in both models (Sect. 3), including the connections with the surface winds and buoyancy forcing, and their transient adjustment. Control of the atmospheric CO 2 by variations in ocean carbon inventories is identified (Sect. 4) and, finally, the wider implications discussed (Sect. 5).

Model design of experiments
To test the skill of our metric, suites of coupled physical and biogeochemical ocean model experiments are performed with parameterized ocean eddies integrated globally (Lauderdale et al. 2013) or with explicit eddies integrated over a sector (Munday et al. 2013(Munday et al. , 2014. It is not yet computationally feasible to conduct global, eddyresolving model integrations to a long-term equilibrium for the carbon system. Both sets of experiments employ the MIT ocean general circulation model  with online biogeochemistry (e.g. Dutkiewicz et al. 2006). The experiments include a well-mixed slab atmosphere, either taking up or returning CO 2 to the ocean, and all the experiments are integrated for 10,000 years and reach a global equilibrium. These experiments, employing different domains and resolutions, are designed to reveal the extent that Southern Ocean dynamics controls atmospheric CO 2 for a long-term equilibrium.

Model experiments with parameterized eddies
A coarse-resolution MITgcm configuration is used for the global ocean with a realistic domain and grid spacing of 2.8° × 2.8° and 15 non-uniform levels in the vertical. Such coarse grid spacing allows the model to be integrated to steady state with different forcing and eddy closure schemes as described below. The model is forced by a 12-monthly repeating cycle of heat and freshwater fluxes and wind stresses with additional relaxation to surface temperature and salinity with a freshwater/salt conserving scheme (Lauderdale et al. 2013). Changes in wind forcing generate associated changes in surface buoyance forcing. The biogeochemistry model carries dissolved inorganic carbon, alkalinity, organic and inorganic phosphate, oxygen and iron tracers (Dutkiewicz et al. 2006); see Lauderdale et al. (2013) for further configuration details.
In the majority of experiments of Lauderdale et al. (2013), a constant eddy thickness diffusivity, K GM , of 1000 m 2 s −1 is employed to parameterize the slumping of density surfaces in baroclinic instability (Gent and McWilliams 1990;Danabasoglu et al. 1994;Gent et al. 1995). The combination of surface wind and buoyancy forcing then leads to a time-mean Eulerian circulation, which is partly opposed by the parameterised eddy circulation, and leads to the residual circulation.
For the ensemble of experiments presented here, several methods are implemented to parameterize the effects of ocean eddies. The eddy thickness diffusivity is allowed to vary horizontally or with depth using three different choices of eddy closure to provide a more realistic representation of eddy exchange. Firstly, K GM is evaluated as proportional to the square of an eddy transfer length scale and the Eady growth rate averaged between depths of 100-2000 m (Visbeck et al. 1997), where the proportionality coefficient α = 0.013 and the eddy transfer length scale l = 200 km, f is the Coriolis parameter and Ri is the Richardson number; these choices are equivalent to setting α = 0.13 and choosing the eddy transfer length scale to be the same as the mid-latitude baroclinic Rossby radius of deformation (l ~ 60 km, Chelton et al. 1998) following Bryan et al. (1999). This thickness diffusivity is also bounded to lie within the range 300-4500 m 2 s −1 . Secondly, the thickness diffusivity is varied in the vertical with the buoyancy frequency (Ferreira et al. 2005), where K ref is set to 4500 m 2 s −1 in the surface layer (e.g. Danabasoglu and Marshall 2007) and N ref 2 is the surface value of the buoyancy frequency N 2 . Thirdly, K GM is allowed to vary in both the horizontal and vertical based upon the stratification and vertical velocity shear from the Richardson number (Hofmann and Morales Maqueda 2011).
Experiments using these different eddy closures for thickness diffusivity are initialized from a previously spunup state and integrated for 10,000 years using the default monthly wind stresses (Lauderdale et al. 2013). From each of these separate control states, four wind stress perturbations (50 % increase in peak Southern Ocean winds, 50 % decrease in peak winds, 10° northward shift of peak winds and 10° southward shift of peak winds) are integrated to steady state for a further 10,000 years (Table 1). Two further experiments with combined perturbations are included, first with peak winds both shifted 3° northward and decreased in strength by 50 % and second with peak winds The surface forcing, summarized by the latitude and peak strength of the Southern Hemisphere westerlies, drives anomalies in atmospheric CO 2 through changes in the wind-driven Eulerian-mean (ψ Eul ) and parameterized eddy (ψ eddy ) circulations that partially compensate to form the Southern Ocean residual circulation (ψ res , maximum value below 500 m). The difference in strength of the upper cell (as measured by ψ res ) and the strength of the circulation associated with the bottom cell (ψ bot , maximum value below 500 m) makes the residual upwelling, as defined by (3). Changes in the circulation affects nutrient distributions, causing changes in the net rate of biological production and the ratio of regenerated nutrient to the total nutrient concentrations evaluated by P*, measuring the efficiency of the soft tissue pump (Ito and Follows 2005)  shifted 3° southward and increased in strength by 50 % (Lauderdale et al. 2012(Lauderdale et al. , 2013. The latter of these perturbations resembles current trends in the Southern Hemisphere westerlies as diagnosed from atmospheric reanalyses (e.g. Huang et al. 2006). Each of these different wind forcings are associated with different buoyancy forcings via the surface relaxation in the model. In all of the model simulations, the total atmosphere-ocean inventory of carbon is conserved, so that only the atmosphere-ocean partitioning of carbon is altered.

Eddying model experiments
The eddying model employs an idealized flat-bottomed sector configuration of MITgcm, 20° wide and 5000 m deep, extending from 60°N to 60°S. The domain includes a re-entrant channel south of 40°S, partially blocked by a half-depth abyssal barrier representing Drake Passage. The model has a 1/2° grid spacing in the horizontal and 42 unevenly spaced levels in the vertical (the surface layer extending over the upper 10 m and the bottom layer over the lower 250 m). The inter-hemispheric basin allows a poleto-pole overturning circulation to be established (Wolfe and Cessi 2011) without relying on equatorial sponge layers or boundary restoring (Munday et al. 2013(Munday et al. , 2014. Surface temperature and salinity are restored to idealized sinusoidal profiles that are elevated at the equator and low at the poles, employing restoring time scales of 10 and 30 days for temperature and salinity respectively. Wind forcing is imposed as an idealized sinusoidal westerly jet between 30°S and 60°S, with the peak at 45°S (off-set slightly north of the channel centre) and zero wind stress north of 30°S. The wind stress peak is varied over the range of 0-1 N m −2 in 12 individual experiments that are each integrated for 10,000 model years; see Munday et al. (2014) for further configuration details. The configuration includes an adiabatic sub-grid scale closure with constant thickness diffusivity, K GM , with a small magnitude of 10 m 2 s −1 (Roberts and Marshall 1998). The majority of tracer transports are carried by the simulated circulation, rather than by parameterized fluxes. For example, the adiabatic parameterisation does not contribute to the net poleward heat transport in the model.
At the northern edge of the re-entrant channel, the deformation radius is ~30-60 km, roughly the same as 1-3 model grid boxes at the horizontal resolution of 1/2°. The model circulation forms closed vortices that are several multiples of the deformation radius. There is an energetic mesoscale circulation, albeit with a smaller eddy kinetic energy than observed via altimetry. This model resolution should be viewed as eddy permitting, as fronts and filaments are not formed that emerge at higher resolution, such as in the 1/6° simulations of Munday et al. (2013). Nevertheless, at 1/2° resolution, circumpolar transport through the re-entrant channel of this sector model is relatively insensitive to wind stress (see Fig. 3a in Munday et al. 2013).
Our goal with this model is to investigate the longterm carbon response to Southern Ocean upwelling due to changes in physical forcing, approaching a thermodynamic and biogeochemical equilibrium, without the response relying upon sub-grid scale eddy closures. In comparison, physical-only integrations to 800 years have been completed at a 1/6° resolution for this sector configuration (Munday et al. 2013) and eddying global ocean model integrations with coupled physics and biogeochemistry have been completed for 200 years at a 1/10° resolution (Dufour et al. 2015).

Model diagnostics
In these model experiments, the residual circulation is defined by a zonal average of the meridional overturning (Marshall and Radko 2003) and includes the partially competing effects of the northward wind-driven Ekman transport and the southward upper ocean transport from geostrophic eddies (Marshall 1997;Gent et al. 1995;Danabasoglu et al. 1994). The residual circulation is the overturning circulation that transports tracers.
The residual upwelling is defined by the total upward transport over the Southern Ocean (displayed in Fig. 1), made up of the difference of the residual circulation across the Southern Ocean and upwelling cells contained south of the ACC that feed the formation of dense bottom waters adjacent to Antarctica. In the control case, the residual upwelling is taken from the difference between the maximum of ψ res between 60°S and 40°S and the minimum of ψ res south of 60°S between 500 and 2000 m depth, defined as The exact latitudinal boundaries are modified in our global model experiments when the westerlies are shifted north or south.
The experiments in both model configurations include a well-mixed slab atmosphere, either taking up or returning CO 2 to the ocean. Both the global and eddying sector model use the same biogeochemical model and all the experiments are integrated for at least 10,000 years and reach a global equilibrium. The global inventory of carbon in the combined atmosphere and ocean is conserved, as the model does not include anthropogenic carbon emissions, riverine carbon input or sediment interactions. (3)

Relation between residual upwelling and atmospheric CO 2
The connection between residual upwelling and atmospheric CO 2 , as depicted schematically in Fig. 1, is now assessed for two different model configurations and with different choices of eddy closures.

Global eddy-parameterized experiments with constant K GM
In our global eddy-parameterized experiments using constant K GM , atmospheric CO 2 increases with stronger westerly winds driving a stronger Eulerian-mean overturning circulation (Fig. 2a), consistent with prior studies (e.g. Lauderdale et al. 2013). However, the eddy-induced circulation might be expected to alter with changing winds: stronger winds steepen potential density surfaces and, in response, a stronger eddy circulation is generated from the flattening of potential density surfaces, and so partly compensate the increase in Eulerian-mean circulation (Viebahn and Eden 2010;Munday et al. 2013). Imposing a doubling of the constant eddy isopycnal thickness diffusivity with stronger winds (e.g. Fyfe et al. 2007) leads to no additional change in the strength of the wind-driven Eulerian-mean circulation, but a slightly reduced increase in atmospheric CO 2 (Fig. 2a, red dots), compared to increased winds with an unchanged value of the isopycnal thickness diffusivity (Fig. 2a, blue dots). Similarly, halving the constant eddy thickness diffusivity with weaker winds does not further alter the strength of the wind-driven Eulerian-mean circulation, but produces a relatively small decrease in atmospheric CO 2 (Fig. 2a, red dots), compared to decreased winds with an unchanged value of the isopycnal thickness diffusivity (Fig. 2a, blue dots). The prescribed increased sensitivity in the eddy-induced circulation acts then to oppose the effects of the wind-driven changes on atmospheric CO 2 . While there is a positive correlation between atmospheric CO 2 and the Eulerian-mean circulation, a tighter relationship is obtained between atmospheric CO 2 and the residual upwelling metric: the slope of the relationship is now relatively insensitive to the change in eddy thickness diffusivity (Fig. 2b). This improved relationship is a consequence of how ocean carbon anomalies are transported by the residual circulation and residual upwelling, rather than by the Eulerian-mean circulation. An increase in the residual upwelling enhances atmospheric CO 2 by reducing the ocean store of carbon. Hence, these experiments with parameterized eddies suggest that the residual upwelling is a useful metric of how Southern Ocean dynamics affects atmospheric CO 2 . This inferred relationship between residual upwelling and atmospheric CO 2 (Fig. 2b) is not reliant on the nature of the eddy closure, illustrated so far only using experiments with different, constant values of thickness diffusivity. The robustness of our metric, and its relationship with equilibrium atmospheric CO 2 , is now demonstrated in two ways by repeating our assessments using global experiments with different, spatially and temporally variable eddy closures and sector eddying experiments, both integrated to equilibrium.

Global eddy-parameterized experiments with variable K GM
Parameterizing mesoscale eddies by employing a non-constant value of the isopycnal thickness diffusion coefficient, K GM , allows coarse resolution models to better resemble the response of eddying models to changing wind forcing, in terms of eddy saturation and more complete eddy compensation (Farneti and Gent 2011;Gent and Danabasoglu 2011;Hofmann and Morales Maqueda 2011). The relationship between residual upwelling and atmospheric CO 2 is now assessed when the eddy thickness diffusivity is allowed to vary, mimicking the likely effect of stronger ocean eddies feeding back on the background ocean state. The eddy diffusivity is chosen either to vary horizontally (constant values in depth) with a measure of the eddy growth rate from baroclinic instability (Visbeck et al. 1997), to change only with depth according to the local buoyancy frequency (Ferreira et al. 2005) or to vary both horizontally and with depth based upon the stratification and vertical velocity shear (Hofmann and Morales Maqueda 2011).
At equilibrium in the control state, westerly winds in the Southern Ocean drive upwelling of deep waters to the surface along potential temperature and density surfaces that strongly shoal towards the south ( Fig. 3 with a peak wind stress of 0.2 N m −2 ). This window between the surface and deep ocean is highlighted by elevated DIC and nutrient concentrations surrounding Antarctica (Figs. 3, 4b, c). Elsewhere, there are relatively low concentrations of DIC and nutrients in the upper ocean due to relatively warm waters and biological utilisation, which become elevated at depth due to the accumulated regeneration of biological fallout. There are relatively low concentrations of DIC at mid-depths in the Atlantic associated with recently ventilated North Atlantic Deep Waters (Fig. 3, right section) and instead higher concentrations of DIC at mid-depths in the North Pacific (Fig. 3, left section) due to limited ventilation and accumulation of regenerated carbon.
In the upwelling regions to the south of the ACC (indicated by streamlines in Fig. 4a), there is elevated surface DIC and phosphate (Fig. 4b, c), and associated outgassing of CO 2 to the atmosphere (Fig. 4d, blue). Conversely, there is oceanic CO 2 uptake further north where Southern Ocean waters are subducted into the ocean interior (Fig. 4d, red), as well as some CO 2 uptake adjacent to Antarctica where upwelled waters cool further.
The Southern Ocean residual circulation and upwelling increases when forced by stronger mid-latitude westerly winds (Table 1). The residual circulation increases between 40°S and 70°S (Fig. 5a, b) associated with strengthening of both the upwelling cells. This strengthening is accompanied by an increase in the values of K GM (Fig. 5c) simulating enhanced eddy circulation. These anomalies are mirrored for the cases with 50 % weaker winds (Table 1).
When the Southern Hemisphere westerlies are shifted 10° south, the Southern Ocean residual circulation and upwelling decrease (Table 1 Surface level and zonally-averaged Pacific (left, nominally along the dotted line at 180°W) and Atlantic (right, nominally at 30°E longitude) sections of DIC (mol C m −3 , coloured blue to yellow) after 10,000 years with control wind stress and overlaid potential temperature surfaces (black contours every 1 °C) effect of the increased Coriolis parameter and reduced zonal circumference at higher latitudes, as well as a decoupling of the latitude of the ACC from the latitude of peak wind stress (Lauderdale et al. 2012). The southward shift of the residual overturning cell is associated with an overturning decrease between 20°S and 60°S, and a slight increase further south (Fig. 5e). A closer alignment between the southward-shifted westerly wind belt and the unblocked latitudes of Drake Passage leads to the overturning return flow of the upper cell being at depth below the sill. In contrast, for northward-shifted winds, the maximum in the overturning streamfunction is far from Drake Passage and the return flow is shallower and intensified by the relatively weak Coriolis parameter and the increased zonal extent of the ocean at lower latitude. A southward shift in westerly winds leads to associated anomalies in K GM (Fig. 5f) with higher values in the vicinity of steepened isopycnals, such as south of 60°S.
The two experiments involving combined changes in westerly wind strength and a small shift in latitude fall within the prior range of model responses. When the wind stress is increased by 50 % and shifted 3° south, the response in the Southern Ocean is a slightly weaker residual upwelling and smaller increase in atmospheric CO 2 , compared with for an increase in wind-stress magnitude alone. Similarly, for 50 % reduced winds with a 3° northward shift, the response is a slightly stronger residual upwelling and elevated atmospheric CO 2 .
The simulated values for the thickness diffusivity compare favourably to those from other studies: Hofmann and Morales Maqueda (2011) doubled Southern Ocean winds in the Modular Ocean Model (MOM3) to obtain increases in K GM of the order 300 m 2 s −1 , as compared with our increase in K GM of ~100-150 m 2 s −1 for 50 % stronger winds (Fig. 5c). For the other parameterizations, Danabasoglu and Marshall (2007) implemented the Ferreira et al. (2005) scheme into the Parallel Ocean Program (POP1.4) obtaining a rapid decay of K GM from 4000 to 400 m 2 s −1 in the upper 500 m near the equator and at high latitudes, while K GM typically varied from 400 to 500 m 2 s −1 at depths of 2000-3000 m between 20°S and 60°S. Finally, Farneti and Gent (2011) added a Visbeck et al. (1997) parameterization to the 1° GFDL climate model and found K GM typically reaching 600 m 2 s −1 between 40°S and 60°S.
The atmospheric CO 2 again systematically increases, nearly linearly, with increasing Southern Ocean residual upwelling (Fig. 6a, b), at an ensemble-averaged rate of 1.22 parts per million/10 6 m 3 s −1 or 1.22 ppm for every 1 Sv (1 Sv ≡ 10 6 m 3 s −1 ; Table 2). There is less scatter in the atmospheric CO 2 relationship with residual upwelling, than with residual overturning, due to the additional accounting of the upwelling south of the ACC. In particular, the model response for decreased and northward-shifted winds includes the effect of changes in the polar easterly winds south of the ACC (Stewart and Thompson 2012).

Eddy-permitting experiments with a sector model
The Southern Ocean residual upwelling metric is diagnosed from a suite of idealized sector model integrations in which mesoscale eddies are simulated explicitly. The model also includes a sub-grid scale adiabatic parameterization, but this closure only has a minor effect on the transport of tracers at this resolution. The eddying integrations of the sector model are integrated for 10,000 years with online biogeochemistry, which is sufficiently long for the model to reach thermodynamic and biogeochemical equilibrium.
At dynamic equilibrium, there is elevated eddy activity near both the northern and southern boundaries ( the sensitivities of CO 2 to residual upwelling are 1.11 ppm Sv −1 (R 2 = 0.97) for the blue, 1.40 ppm Sv −1 (R 2 = 0.94) for the orange and 1.15 ppm Sv −1 (R 2 = 0.98) for the red points. For the ocean sector model integrations with eddying circulation, the red in (c), (d), the sensitivity of CO 2 to residual upwelling is 8.46 ppm Sv −1 (R 2 = 0.98) for a peak wind stress of 0.2 N m −2 ), where potential temperature and density surfaces strongly shoal poleward. The surface wind stress generates more eddy kinetic energy in the re-entrant channel. There is a vigorous meridional overturning cell providing well ventilated waters with relatively low DIC in the upper waters above the sill in the re-entrant channel. The lack of wind forcing north of 30°S prevents gyre circulations forming, so that isotherms are flat in the low latitudes of the closed basin. The highest concentrations of DIC are in the abyss near the northern boundary due to limited ventilation.
Eddy signatures are revealed in the horizontal transport streamfunction (Fig. 8a) and undulations of tracer distributions, such as for potential temperature (Fig. 7), DIC and phosphate (Fig. 8b, c), particularly evident at high latitudes. There is also an imprint of eddy activity on the air-sea flux of carbon (Fig. 8d), despite the air-sea equilibration timescale being typically a year.
When the wind stress is increased in the eddying model, the increased wind drives an increase in eddy kinetic energy and results in eddy saturation of the circumpolar transport and partial eddy compensation of the overturning   (Munday et al. 2013). As with the global model, atmospheric CO 2 again increases with stronger residual circulation and residual upwelling (Fig. 6c, d). Again, in common with the global model, there is a tighter relationship with a more linear response for the atmospheric CO 2 versus residual upwelling (Fig. 6d): an increase in residual upwelling leading to higher atmospheric CO 2 due to enhanced communication between the deep ocean and atmosphere. The exception to this linear response is close to where the residual upwelling becomes vanishingly small and there is very low atmospheric CO 2 ; here the deep ocean becomes very isolated and the carbon balances moves to a diffusive regime.

Connection to surface buoyancy forcing in the global ocean experiments
Our model experiments automatically include changes in both wind and buoyancy forcing (Fig. 9a, b), the latter involving changes in air-sea heat and freshwater fluxes via a surface relaxation in the model. Changes in wind forcing alter the northward Eulerian-mean circulation across the Southern Ocean, part of which is compensated by mesoscale eddies, and thus alter the residual circulation and upwelling. From the suite of global model experiments, atmospheric CO 2 is generally positively correlated with the zonal wind strength (Fig. 9a), but the slope is dependent on the strength of the feedback from the eddy circulation; hence, there is no simple relationship between these variables. The residual circulation and residual upwelling only occur if there is an associated buoyancy forcing driving a transformation of dense waters. There is a generally linear relationship between atmospheric CO 2 and the buoyancy forcing integrated over ACC streamlines (Fig. 9b), but there is also significant scatter using this simple measure of buoyancy forcing. The residual upwelling may also be defined from the convergence in the transformation of water masses controlled by buoyancy forcing (Fig. 9c), based upon the divergence in water-mass transformation rate, G, or diapycnal volume transport, where σ is the potential density minus 1000 kg m −3 . The water-mass transformation rate, G, is derived from the buoyancy flux across potential density outcrops marking the northern flank of the ACC and off the Antarctic shelf break (Walin 1982;Marshall 1997;Speer et al. 2000;Badin and Williams 2010;Badin et al. 2013). For the control and the model experiments with changing strength in the winds, this convergence is evaluated as the peak transformation to denser waters at σ = 27.4 minus the transformation to lighter surfaces at σ = 27.0 using (4); different bounding potential density surfaces are chosen according to their geometry when the westerlies are shifted in latitude. The atmospheric CO 2 again systematically increases with the Southern Ocean buoyancy-derived residual upwelling, at a rate of typically 1.61 parts per million/10 6 m 3 s −1 or 1.61 ppm for every 1 Sv (Table 2). This dependence between the atmospheric CO 2 and the buoyancy-derived residual upwelling is similar to that derived for the kinematic diagnostics of the residual upwelling (Fig. 6b). The slight differences are due to differences in upwelling estimates. The buoyancy-derived estimate of residual upwelling should also account for the effects of diapycnal mixing in water-mass transformation. However, the diapycnal mixing contribution to transformation in the Southern Ocean is diagnosed to make a minor contribution for lighter mode waters, but becomes more significant for denser deep and bottom waters (Badin et al. 2013).

Transient adjustment of atmospheric CO 2 relationship in the global ocean experiments
There is a transient adjustment in our experiments, initially with a weak positive relationship between atmospheric CO 2 and residual upwelling, which steepens and becomes less scattered until an equilibrium slope is approached after 500 years (Fig. 10). Hence, our metric of residual upwelling is more important in determining atmospheric CO 2 on timescales from several centuries to millennia. For our suite of experiments, the equilibrium slope is close to an atmospheric CO 2 increase of 1.5 ppm for every 1 Sv increase in residual upwelling. (x10 Fig. 9 The effect of the Southern Ocean on atmospheric CO 2 based on global model experiments with parameterised eddies after 10,000 years. Scatterplots for atmospheric CO 2 mixing ratio (ppm) versus a the root-mean square zonal wind stress (N m −2 ), b the buoyancy forcing over streamlines in the Antarctic Circumpolar Current (10 −9 m −2 s −3 ) and c the buoyancy-derived upwelling strength (Sv).
The model experiments range from weak to strong residual circulation and upwelling and employ 3 different eddy closures with variable K GM from Visbeck et al. (1997)  A dependence of atmospheric CO 2 to residual upwelling of 1 ppm versus 1 Sv change is denoted by the dashed line. The sensitivities of CO 2 to residual upwelling range for the different eddy closures from 0.11-0.15 ppm Sv −1 at 10 years, increasing to 0.51-0.57 ppm Sv −1 at 100 years, 0.79-0.94 ppm Sv −1 at 500 years, 0.88-1.09 ppm Sv −1 at 1000 years, 1.11-1.39 ppm Sv −1 at 5000 years, and 1.11-1.40 ppm Sv −1 after 10,000 years changes in the subduction of mode waters from the winter mixed layer into the thermocline. On centennial and longer timescales, the global thermocline is also likely to adjust through altered Southern Ocean overturning (e.g. Allison et al. 2011;Jones et al. 2011).

Relationship between atmospheric CO 2 changes and changes in ocean carbon pools in the global ocean experiments
The connection between atmospheric CO 2 and the ocean carbon cycle is now explored by diagnosing the different ocean carbon reservoirs. The ocean carbon changes are understood by separating the DIC into different carbon pools (Ito and Follows 2005;Williams and Follows 2011;Lauderdale et al. 2013;Munday et al. 2014), where C sat is the saturated concentration of DIC at equilibrium, C soft and C carb are the concentration of DIC resulting from regeneration of biogenic soft-tissue and carbonate and C res is the disequilibrium concentration, which is the residual between DIC and C sat when the water mass was last at the sea surface. Changes in C sat may be further separated into four components to reveal the change in concentration of carbon due to each contributing factor (6), assuming a linearization of the buffering capacity of the ocean carbonate system (Goodwin and Lenton 2009), which is a reasonable assumption over the range of atmospheric CO 2 values considered here.
where the subscripts on the differentials represent those variables being kept constant. It is convenient to consider the forced effects on the saturation carbon reservoir caused by changes in potential temperature (θ), salinity (S) and preformed alkalinity (A pre ) separately from the negative feedback due to how the saturated carbon alters with atmospheric pCO 2 . Furthermore, due to the strongly compensating interactions between C soft and C res highlighted by Lauderdale et al. (2013), these reservoirs are combined together as the "non-saturated" component, C nonsat , of DIC, giving a resulting three-component separation, In our experiments, atmospheric CO 2 changes are mirrored by compensating changes in the ocean carbon DIC = C sat (θ, S, A pre ) + C sat (pCO 2 ) + C nonsat .
reservoirs (Ito and Follows 2005;Williams and Follows 2011;Lauderdale et al. 2013;Munday et al. 2014). For example, an increase in the residual upwelling from stronger Southern Ocean winds is associated with DIC concentrations increasing in the upper ocean and in deep waters at high latitudes, but decreasing elsewhere over much of the ocean (Fig. 11a). This response is made up of contrasting contributions from the different carbon pools. Firstly, at equilibrium under a higher atmospheric CO 2 , the surface ocean carbon anomalies increase to enable the globally-averaged air-sea gas exchange to vanish. This response is partly achieved by a strong negative feedback with the ocean saturated carbon anomaly associated with pCO 2 increasing with rising atmospheric CO 2 (Fig. 11b). Secondly, stronger winds and increased residual upwelling consistently decrease the ocean saturated carbon pool associated with temperature, alkalinity and salinity changes (Fig. 11c), which is to be expected from a thicker and warmer thermocline with enhanced subduction of intermediate and mode waters (Lauderdale et al. 2013). Thirdly, this decrease in the saturated carbon pool is reinforced by a reduction in the non-saturated carbon pool (Fig. 11d). Stronger residual upwelling and enhanced subduction increases the proportion of preformed nutrients in the ocean interior and, thus, reduces the ocean carbon store by decreasing the efficiency of the biological pump; as indicated by generally low P* values (Table 1) where P* is the average ratio of regenerated nutrient to total nutrient concentrations (Ito and Follows 2005). This decrease in the ocean carbon store occurs despite an increase in global biological production driven by upwelled nutrients (Lauderdale et al. 2013). The non-saturated carbon anomalies involve large partially opposing anomalies in the regenerated carbon and the disequilibrium carbon: increased upwelling of regenerated carbon and DIC from the deep Southern Ocean is incompletely outgassed to the atmosphere with a reduced surface residence time (Lauderdale et al. 2013).
The changes in ocean carbon components (Fig. 12) generally follow a similar response across our model experiments with the three sets of eddy closures: increased winds and/or northward-shifted winds lead to enhanced residual upwelling and the carbon components varying to provide elevated atmospheric CO 2 , while decreased winds lead to decreased residual upwelling and reduced atmospheric CO 2 .
However, there is a different response when the winds are shifted southward, leading to weaker residual upwelling and the change in carbon components leading to lower atmospheric CO 2 . The DIC anomaly (Fig. 11e) in intermediate depths is increased both through a large positive anomaly in the saturated carbon pool from temperature, alkalinity and salinity (Fig. 11g), due to a shoaled , which is the sum of opposing biogenic regenerated carbon and disequilibrium carbon concentrations pycnocline and relatively cooler upper ocean, and a positive anomaly in the upper ocean non-saturated carbon pool (Fig. 11h), driven by an increase in the biologically-regenerated carbon. Even though the rate of residual upwelling is lower (Table 1), the southward-shift in winds drives an upwelling of denser isopycnals with more regenerated carbon, which leads to a negative anomaly in the deep ocean C soft pool and low values for P* (Table 1). At the surface these carbon anomalies are only partially dampened by airsea outgassing before the remaining regenerated carbon is subducted into the ocean interior as a positive disequilibrium carbon anomaly. The overall effect is that increasing residual upwelling leads to a rise in atmospheric CO 2 by reducing the ocean store of carbon, which is achieved by reducing the saturated carbon pool associated with changes in temperature, alkalinity and salinity, and the non-saturated carbon pool (Fig. 12).
There are similar carbon inventory changes in the eddying sector model integrations (Fig. 12d). When surface wind stress is increased, there is an increase in the atmospheric CO 2 due to a combined decrease in both the saturated and non-saturated ocean carbon pools, in common with the global model. When surface wind stress is decreased, there is a slight decrease in the atmospheric CO 2 , but there are partly compensating changes in both the saturated and non-saturated ocean carbon pools. The detailed response of the individual ocean carbon pools may differ then with model resolution and complexity, even though their integrated response may be similar.

Discussion
This modelling study suggests that residual upwelling provides a robust metric of how ocean dynamics in the Southern Ocean affect atmospheric CO 2 , giving an increase in atmospheric CO 2 of typically 1.5 ppm for every 1 Sv increase of residual upwelling in the Southern Ocean.

Implications for the Southern Ocean control of atmospheric CO 2
Atmospheric CO 2 changes have been separately explained in terms of Southern Ocean changes in wind (Toggweiler et al. 2006;Le Quéré et al. 2007;Tschumi et al. 2008;Lovenduski et al. 2013Landschützer et al. 2015) and buoyancy forcing (Marshall 1997;Watson and Naveira Garabato 2006;Watson et al. 2015). Our metric of residual upwelling provides a unified view by accommodating the effects of both these forcings. The residual circulation across the Southern Ocean is directly linked to buoyancy forcing (Marshall 1997;Speer et al. 2000;Badin and Williams 2010;Morrison et al. 2011) and the formation of water masses (Walin, 1982;Speer et al. 2000;Badin et al. 2013). The wind and buoyancy forcing are tightly coupled, as the buoyancy forcing may adjust to the wind and eddydriven transfer across the ACC (Badin and Williams 2010). The suite of experiments with parameterized and explicit eddies both support the view that residual upwelling is an important metric of how the Southern Ocean controls atmospheric CO 2 on centennial to millennial timescales. Previous model studies have suggested that enhanced communication between the surface and deep Southern Ocean leads to elevated atmospheric CO 2 , either employing box models (Knox and McElroy 1984;Sarmiento and Toggweiler 1984;Siegenthaler and Wenk 1984;Toggweiler 1999) or more dynamically-complete circulation models of the Southern Ocean (Toggweiler et al. 2006;Lovenduski et al. 2007Lovenduski et al. , 2008Menviel et al. 2008;Tschumi et al. 2008Tschumi et al. , 2011Lauderdale et al. 2013;Munday et al. 2014). These prior studies have emphasized either the effect of the Eulerian-mean upwelling or the residual overturning circulation. In contrast, our metric explicitly accounts for the effect of ocean overturning changes in both the upper and lower cells, including partially-compensating effects of the parameterized eddy circulation.
The sensitivity of atmospheric CO 2 to residual upwelling is much larger in the sector model than the global model. This reduced sensitivity in the global model is probably due to the redistribution of carbon between additional basins, such as the Pacific, which reduces the amount of carbon that may be exchanged with the atmosphere. These two models of differing complexity employed here represent different limits of the spread of climate models. There is a broadly similar climate response in the Southern Ocean to forcing in the Coupled Model Inter-comparison Project phase 3 (CMIP3) and the more advanced CMIP5 integrations (Meijers 2014). Hence, we anticipate that this linear relationship between residual upwelling and atmospheric CO 2 will carry over to more complex climate models and this relationship might be viewed an emergent property of climate models.
In summary, increasing residual upwelling leads to a rise in atmospheric CO 2 by increasing the communication with the deep ocean and reducing the ocean store of carbon. This linear response is found in two ocean models with different domains and representations of ocean eddies. Within the global model calculations, the long-term effect of Southern Ocean dynamics is that an increase in residual upwelling of 1 Sv corresponds to an increase in atmospheric CO 2 of typically 1-1.5 ppm. In principle, the residual upwelling might be diagnosed through repeated estimates of the meridional velocity and density in the Southern Ocean to calculate (3) or by evaluating the watermass transformation (4) from estimates of the surface density, air-sea heat and freshwater fluxes, and preferably together with estimates of diapycnic mixing.

Implications for the paleoceanographic context
Enhanced Southern Ocean upwelling is inferred from greater opal accumulation due to enhanced supply of dissolved silica during the transition from glacial to interglacial climate (Anderson et al. 2009). This tracer-based inference of enhanced upwelling is consistent with our residual upwelling metric, as residual upwelling encapsulates both the upper and lower circulation cells (e.g. in Fig. 1), rather than focusing on only one of those cells involving either Eulerian-mean upwelling or residual circulation.  Fig. 12 The effect of the Southern Ocean dynamics on the ocean carbon pools (left axis) ordered by change in atmospheric CO 2 (50 % stronger westerly winds, 10° northward-shifted winds, 50 % stronger westerly winds shifted 3° southward, 10° southward-shifted winds, 50 % weaker winds shifted 3° northwards, and 50 % weaker magnitude winds), which correlates with the change in residual upwelling (diamonds, right axis). The global circulation ocean model integrations for 10,000 years have parameterised eddies using 3 different eddy closures of variable K GM from a Visbeck et al. (1997), b Ferreira et al. (2005 and c Hofmann and Morales Maqueda (2011). Car-bon partitioning from the eddy-permitting sector model is shown in (d). The change in atmospheric CO 2 anomaly (dark blue bar) results from changes in the ocean carbon inventory, which is separated into a saturated carbon pool from CO 2 partial pressure (pCO 2 ) changes (light blue bar), non-pCO 2 changes (from temperature, salinity and alkalinity, orange bar) and non-saturated carbon changes (red bar) from regenerated carbon and disequilibrium carbon. The increase or decrease in atmospheric CO 2 is generally achieved via opposing increases or decreases in saturated carbon from non-pCO 2 changes (orange bar) and/or non-saturated carbon (red bar) The atmospheric CO 2 change between glacial and interglacial climates is typically 80-100 ppm (Petit et al. 1999). Our diagnostics of residual upwelling suggest an increase in atmospheric CO 2 of typically 1.5 ppm for every 1 Sv change in residual upwelling. An upper bound of how important Southern Ocean upwelling is for atmospheric CO 2 is then provided by the maximum possible change in the residual upwelling multiplied by our sensitivity in atmospheric CO 2 . Taking the maximum possible circulation change to be the modern overturning of 20-25 Sv and collapsing to 0 Sv, then suggests a maximum change in atmospheric CO 2 of typically 30-40 ppm.
This upper bound estimate does not make any assumptions about how the change in residual upwelling is accomplished, because the same sensitivity could be achieved by winds or buoyancy forcing ( Table 2). The idealized experiments presented here have merely used Southern Ocean winds as a convenient way to alter the residual upwelling. Indeed, westerly wind changes during the last glacial-interglacial transition were probably more modest than used in our experiments (Kohfeld et al. 2013;Sime et al. 2013) and so buoyancy fluxes could provide supplementary physical forcing (e.g. Watson and Naveira Garabato 2006;Ferrari et al. 2014;Watson et al. 2015).
Our upper bound of 30-40 ppm is comparable with the estimate of a 27 ppm contribution to the glacial-interglacial CO 2 change from "altered ocean circulation" albeit within wide uncertainty bounds of 3-57 ppm (Kohfeld and Ridgwell 2009). Hence, other processes need to be invoked to explain the rest of the 80-100 ppm CO 2 anomaly over a glacial-interglacial transition, such as changes in biological activity (e.g. Ito and Follows 2013) that might enhance ocean regenerated carbon storage.
Together with other metrics, such as changes in seaice cover (Stephens and Keeling 2000;Ferrari et al. 2014) and the proportion of preformed nutrients (Ito and Follows 2005), our metric of Southern Ocean residual upwelling dependence may be useful in providing bounds as to the likely climatic effect of the Southern Ocean circulation for either the present day or glacial-interglacial cycles.