Marine carbon cycle response to a warmer Southern Ocean: the case of the last interglacial

Recent studies investigating future warming scenarios have shown that the ocean carbon sink will weaken over the coming century due to ocean warming and changes in oceanic circulation. However, significant uncertainties remain regarding the magnitude of the oceanic carbon cycle response to warming. Here, we investigate the Southern Ocean's (SO, south of 40^∘ S) carbon cycle response to warmer conditions, as simulated under last interglacial boundary conditions (LIG, 129–115 ka). We find a ∼150 % increase in carbon dioxide (CO₂) outgassing over the SO at the LIG compared to pre-industrial conditions (PI), due to a 0.5 ^∘C increase in SO sea surface temperatures. This is partly compensated for by an equatorward shift of the Southern Hemisphere (SH) westerlies and weaker Antarctic Bottom Water formation, which both lead to an increase in dissolved inorganic carbon (DIC) in the deep ocean at the LIG compared to PI. These deep-ocean DIC changes arise from increased deep- and bottom-water residence times and higher remineralization rates due to higher temperatures. While our LIG simulation features a large reduction in SO sea ice compared to the PI, we find that changes in sea ice extent exert a minor control on the marine carbon cycle. The projected future strengthening and poleward shift of the SH westerlies coupled to warmer conditions at the surface of the SO should thus weaken the capacity of the SO to absorb anthropogenic CO₂ over the coming century.

1 Introduction

Future increases in atmospheric carbon dioxide (CO₂) concentration are unequivocally projected to further warm the Southern Ocean (SO) and reduce sea ice concentrations (Bracegirdle et al., 2020). The current state of knowledge suggests the mitigating effects of carbon cycle feedbacks on global warming to be less efficient under future scenarios (Pachauri et al., 2014). This primarily stems from the changes in carbon uptake by the terrestrial biosphere and ocean under a changing climate. Both land and ocean presently act as sinks of anthropogenic carbon, each absorbing about 25 % of anthropogenic emissions (Friedlingstein et al., 2020), with 40 % of the ocean sink being attributed to the SO (Caldeira and Duffy, 2000; DeVries, 2014). The SO CO₂ uptake weakened during the 1990s due to a strengthening and poleward shift of the Southern Hemisphere (SH) mid-latitude westerlies (Lovenduski et al., 2007; Le Quéré et al., 2007; Zickfeld et al., 2007; Gruber et al., 2019a, b) but strengthened in the 2000s due to cooling over the Pacific and increased stratification over the Atlantic and Indian sectors of the SO (Landschuetzer et al., 2016; Gruber et al., 2019b).

In addition, the consensus amongst studies analysing future climate simulations points towards amplified warming over high latitudes (Holland and Bitz, 2003; Smith et al., 2019; Fan et al., 2020). For example, the SO annual mean sea surface temperature (SST) anomaly is projected to exceed 0.5 ^∘C for Shared Socioeconomic Pathway scenario (SSP) 245 and 1.5 ^∘C for SSP 585 at 2100 relative to 2015 (Bracegirdle et al., 2020; Tonelli et al., 2021), with reduced sea ice during austral spring (Roach et al., 2020). At the same time, SH westerlies are projected to strengthen and shift poleward over the coming century (Collins et al., 2013; Zheng et al., 2013; Goyal et al., 2021). Overall, climate change will reduce the amount of anthropogenic carbon taken up by the ocean (Plattner et al., 2001; Bernardello et al., 2014; Wang et al., 2014; Howes et al., 2015) by reducing CO₂ solubility (Bernardello et al., 2014), weakening the efficiency of the biological pump (Boyd, 2015), increasing outgassing from upwelling of Circumpolar Deep Waters (CDW) due to stronger and poleward-shifted SH westerlies (Lovenduski et al., 2007; Zheng et al., 2013; Gruber et al., 2019a), and reducing sea ice extent (Shadwick et al., 2021). However, uncertainties still remain in regards to the response of the SO carbon cycle under a warmer climate.

The last interglacial (LIG, 129–115 ka) was the warmest interglacial of the last 800 000 years (Masson-Delmotte et al., 2013; PAGES, 2016). The warmer climate at the LIG is primarily attributed to a stronger Northern Hemisphere summer insolation (Laskar et al., 2004) owing to the orbital configuration of higher eccentricity and obliquity (Berger, 1978), rather than higher greenhouse gas concentrations as projected for the future. The role of orbital forcing versus greenhouse gases on temperature has been analysed in detail in Yin and Berger (2012). The LIG is associated with annual mean SSTs around 0.5 ^∘C higher than pre-industrial conditions (PI) (Capron et al., 2014; Hoffman et al., 2017), and sea level was ∼ 5 m higher than PI, with estimates ranging from 1.2 to 9 m higher (Dutton et al., 2015; Rohling et al., 2019; Dyer et al., 2021). The area-weighted summer warming is estimated to be 1.1–1.9 ^∘C for the North Atlantic and 1.6–1.8 ^∘C for the SO compared to PI (Govin et al., 2012; Capron et al., 2014; Hoffman et al., 2017). Summers over land areas are reconstructed to be 4–5 ^∘C warmer at high latitudes in the Northern Hemisphere (CAPE-Last Interglacial Project Members, 2006), 3–11 ^∘C over Greenland (Dahl-Jensen et al., 2013), and 2.2 ^∘C over Antarctica (Masson-Delmotte et al., 2013). Reconstructions of sea ice are mostly focused on the Arctic and suggest ice-free conditions south of 78^∘ N (Van Nieuwenhove et al., 2008, 2011, 2013; Kageyama et al., 2021). The Atlantic Meridional Overturning Circulation (AMOC) has been suggested to have weakened at the peak of the LIG (127 ka) and strengthened afterwards (Galaasen et al., 2014), with strong evidence pointing towards periods of reduced Antarctic Bottom Water (AABW) formation during the early part of the LIG due to discharges from the Antarctic ice sheet (Hayes et al., 2014; Rohling et al., 2019).

A few studies have investigated the terrestrial carbon response to LIG conditions (Kleinen et al., 2015; Brovkin et al., 2016), but the marine carbon cycle response at the LIG, and particularly in the SO, has received little attention. A comprehensive study of the changes in the SO carbon cycle at the LIG compared to PI can enhance our understanding of the processes involved in the SO carbon cycle and their sensitivity to changes in boundary conditions. It can also help us, to a certain extent, to better quantify the impact of expected future physical and dynamical changes in a warming SO on the marine carbon cycle, bearing in mind that there was no additional anthropogenic carbon released during the LIG.

Here, we study the impact of a warmer climate, in particular at high latitudes, on the oceanic carbon cycle by analysing an equilibrium last interglacial simulation (lig127k, Otto-Bliesner et al., 2017) performed with the ACCESS-ESM1.5 model (Ziehn et al., 2020; Yeung et al., 2021). The paper is structured as follows. In Sect. 2, components of the ACCESS-ESM1.5 model are described, followed by a framework for decomposing variables relevant for carbon cycle changes. The main results are presented in Sect. 3, discussing the changes in climate and the carbon cycle. These include differences in ocean properties, such as temperature, sea ice cover, stratification and circulation, and their effects on the carbon cycle, including air–sea fluxes, net carbon storage, and the efficiency of the biological pump. Finally, in Sect. 4, we conclude and discuss the limitations of the current study, including possible sources of uncertainty, and what inferences can be drawn for future global warming scenarios.

2 Methods

2.1 Model description and experimental design

An equilibrium last interglacial simulation (lig127k) is performed with the Australian Community Climate and Earth System Simulator Earth System Model, ACCESS-ESM1.5 (Ziehn et al., 2020), which includes interactive land and ocean carbon cycles and is Australia's submission to the Paleoclimate Modeling Intercomparison Project 4 (PMIP4) – Coupled Model Intercomparison Project (CMIP6) (Yeung et al., 2021). ACCESS-ESM1.5 differs from the previous version, ACCESS-ESM1 (Law et al., 2017) mostly in the land and ocean components. ACCESS-ESM1.5 is built upon the ACCESS1.4 physical model and includes the UK Met Office Unified Model (UM) version 7.3 (Martin et al., 2010; The HadGEM2 Development Team, 2011) as an atmospheric component, which is directly coupled to the updated land Community Atmosphere Biosphere Land Exchange model (CABLE) version 2.4 (Kowalczyk et al., 2013), both with a horizontal resolution of 1.875^∘ × 1.25^∘. UM has 38 vertical levels. The ocean model is the NOAA/GFDL Modular Ocean Model (MOM) version 5 (Griffies, 2014), which is coupled to the Los Alamos National Laboratory sea ice model (LANL CICE) version 4.1 (Hunke et al., 2010). MOM5's resolution is 1^∘ × 1^∘ with 50 vertical levels. The coupler is the Ocean Atmosphere Sea Ice Soil – Model Coupling Toolkit (OASIS-MCT) (Craig et al., 2017). The current version of the land model CABLE also includes biogeochemistry (BGC) implemented using the CASA-CNP module (Wang et al., 2010). The ocean carbon cycle is simulated using the Whole Ocean Model of Biogeochemistry And Trophic-dynamics (WOMBAT) model (Oke et al., 2013).

WOMBAT is a nutrient–phytoplankton–zooplankton–detritus (NPZD) model, with one class of phytoplankton and one class of zooplankton. It includes dissolved inorganic carbon (DIC), alkalinity (ALK), phosphate (PO₄), oxygen (O₂), and iron. The biogeochemical tracers are coupled using the stoichiometric C (carbon):N (nitrogen):P (phosphorous):O₂ ratio of $\mathrm{106}:\mathrm{16}:\mathrm{1}:\mathrm{172}$. The air–sea gas exchange is based on the square of the wind speed (Wanninkhof, 1992), and the seawater partial pressure of CO₂ (pCO₂) is calculated following the third phase of the Ocean Carbon-Cycle Model Intercomparison Project (OCMIP) protocol using temperature, salinity, DIC, ALK, and PO₄. WOMBAT simulates production in and export from the photic zone and remineralization and dissolution at depth for both organic and inorganic (CaCO₃) particulate matter, with parameters adjusted for inorganic export to be ∼8 % of organic export. The remineralization of organic matter is calculated based on depth- and temperature-dependent parameters, while the dissolution of CaCO₃ uses a constant dissolution rate. ACCESS-ESM1.5 does not include a sediment component, and all particulates reaching the bottom ocean layer are instantly remineralized following the relevant remineralization rates. Other specific details of the model can be found in Oke et al. (2013), Law et al. (2017), and Ziehn et al. (2020).

A pre-industrial 1850 simulation (piControl) is run in accordance with the CMIP6 protocol (Eyring et al., 2016), with a constant atmospheric CO₂ forcing of 284.3 ppm but using CMIP5 solar irradiance (1365.65 Wm⁻², as explained in Ziehn et al., 2020). This simulation is run for 1000 years. Initialized from this PI control run, a LIG simulation (lig127k) is performed using orbital parameters following the PMIP4 protocol (Otto-Bliesner et al., 2017; Yeung et al., 2021) but with the solar constant adjusted to the CMIP5–PMIP3 value to be comparable to the piControl simulation (1365.65 Wm⁻²). CO₂ concentration for the LIG is set at 275 ppm. Vegetation is kept constant at 1850 PI levels. The lig127k experiment is run for 650 years. The run is at equilibrium over the SO, although there is a small drift in globally averaged DIC of ∼1 µmol kg⁻¹/100 yr⁻¹ below 3 km. The analysis presented here is based on the average of the last 100 years of the lig127k simulation compared to the last 100 years of the piControl simulation. Further details on these experiments can be found in Otto-Bliesner et al. (2017), Ziehn et al. (2020), and Yeung et al. (2021).

WOMBAT includes two DIC tracers, one being forced by the prescribed atmospheric CO₂ concentration (PI: 284.3 ppm; LIG: 275 ppm), and the other forced by a constant atmospheric CO₂ concentration of 280 ppm. Unless otherwise stated, all of the analyses presented here are based on the tracers forced with CO₂ concentrations of 280 ppm, while the climate response is forced with the radiative forcing of 284.3 ppm for PI and 275 ppm for the LIG. This allows quantification of the effects of the LIG climate on the carbon cycle independently of the difference in background CO₂ concentrations.

2.1.1 Model evaluation

After a 3000-year and 1000-year spin up of the physical and biogeochemical states, a 500-year piControl run was generated and used to assess the performances of the ACCESS-ESM1.5. The drifts in the piControl simulation are described in Ziehn et al. (2020) and are $-\mathrm{8.5}\times {\mathrm{10}}^{-\mathrm{5}}$ ^∘C century⁻¹ for SST, 5.3 × 10^−3∘C per century for global ocean temperatures, 7.6 × 10⁻⁴ psu per century for sea surface salinity (SSS), and 8.8 × 10⁻⁴ psu per century for global ocean salinity. The drift in total ocean productivity is −0.0163 PgC yr⁻¹ per century and that in carbon flux is −0.0048 PgC yr⁻¹ per century. The net pre-industrial carbon flux is 0.02 and 0.08 PgC yr⁻¹ for land and ocean, respectively.

Ziehn et al. (2020) present results from the historical simulation performed with the ACCESS-ESM1.5 and evaluate the performance with respect to available observations. Here, we summarize the model performance pertaining to the SO. The model has a warm SST bias in the SO, possibly resulting from a too shallow and warm summer mixed layer; however, the climatological sea ice extent in the SO closely follows the observations (Fetterer et al., 2017). The ACCESS-ESM1.5 captures the surface patterns of nutrients relatively well, showing a 0.89 correlation with the surface phosphate distribution from World Ocean Atlas (Garcia et al., 2010b). The simulated primary productivity is higher over 40–45 ^∘S compared to observations (Behrenfeld and Falkowski, 1997), which leads to an increased uptake of CO₂ compared to the observed estimates (Takahashi et al., 2009). Compared to the GLODAP dataset (Key et al., 2004), phosphate and alkalinity concentrations are underestimated in the Southern Ocean by ∼0.4–0.8 and ∼50–100 mmol m⁻³, respectively. In addition, the oxygen concentration is ∼50–100 mmol m⁻³ higher in the ACCESS-ESM1.5 across all depths (Garcia et al., 2010a). These biases could result from a strong ocean ventilation or a weak biological pump in the SO. Further details on model performance and evaluation can be found in Ziehn et al. (2020).

2.2 pCO₂ decomposition

Changes in surface pCO₂, which ultimately control the direction and magnitude of air–sea fluxes, can further be decomposed into contributions from SST, SSS, DIC, and ALK (Sarmiento and Gruber, 2006). The SST contribution is calculated as

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }\mathrm{p}{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{SST}}=e^{({\mathit{\gamma }}_{\mathrm{SST}}\times \mathrm{\Delta }\mathrm{SST})}\times \mathrm{p}{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{PI}}-\mathrm{p}{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{PI}},\end{array}$$

where $\mathrm{\Delta }\mathrm{p}{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{SST}}$ represents the contribution of SST change to surface pCO₂ change and γ_SST=0.0423 ^∘C⁻¹ is the sensitivity of pCO₂ changes to changes in SST (Sarmiento and Gruber, 2006). ΔSST is the difference in SST between LIG and PI, and pCO₂^PI is the surface seawater pCO₂ from the piControl simulation.

The SSS, DIC, and ALK contributions are calculated as

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }{{\mathrm{pCO}}_{\mathrm{2}}}^X={\mathit{\gamma }}_X\times {\displaystyle \frac{\mathrm{\Delta }X}{\stackrel{\mathrm{‾}}{X}}}\times {{\mathrm{pCO}}_{\mathrm{2}}}^{\mathrm{PI}},\end{array}$$

where $\mathrm{\Delta }{{\mathrm{pCO}}_{\mathrm{2}}}^X$ represents the contribution of change in variable X (SSS, DIC, and ALK) to surface pCO₂ change, ΔX is the change in X between LIG and PI, $\stackrel{\mathrm{‾}}{X}$ is the mean value of X at PI, and γ_X is the sensitivity of pCO₂ changes to changes in variable X (γ_DIC is known as the Revelle factor), with γ_SSS=1,

$${\mathit{\gamma }}_{\mathrm{DIC}}=\left\{\begin{array}{lr}\mathrm{15},& \text{if}\mathrm{lat}\ge \mathrm{60}{}^{\circ }\\ \mathrm{13},& \mathrm{lat}<\mathrm{60}{}^{\circ },\end{array}\right.$$

and

$${\mathit{\gamma }}_{\mathrm{ALK}}=\left\{\begin{array}{lr}-\mathrm{12},& \text{if}\mathrm{lat}\ge \mathrm{60}{}^{\circ }\\ -\mathrm{11},& \mathrm{lat}<\mathrm{60}{}^{\circ }.\end{array}\right.$$

These values are based on previous estimates of the meridional profiles of DIC and ALK buffer factors (Sarmiento and Gruber, 2006; Smith and Marotzke, 2008; Egleston et al., 2010; Jiang et al., 2019), although uncertainties still remain regarding these estimates. This decomposition helps shed light on the independent effects of different physical variables to the net surface pCO₂ change.

2.3 Carbon partitioning

To better quantify changes in the carbon cycle between the LIG and PI experiments, we split the total DIC into its remineralized (C_org), dissolved carbonate (C${}_{{\mathrm{CaCO}}_{\mathrm{3}}}$), and preformed (C_pre) components using the equations listed below.

C_org can be estimated from regenerated phosphate (PO${}_{\mathrm{4}}^{\mathrm{Reg}}$) using the C:P stoichiometric ratio (r_C : P=106) (Ito and Follows, 2005):

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }{\mathrm{C}}_{\mathrm{org}}=r_{\mathrm{C}:\mathrm{P}}\times \mathrm{\Delta }{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{Reg}}.\end{array}$$

PO${}_{\mathrm{4}}^{\mathrm{Reg}}$ can be approximated using apparent oxygen utilization (AOU) and the P:O₂ stoichiometric ratio ($r_{\mathrm{P}:{\mathrm{O}}_{\mathrm{2}}}=\mathrm{1}/\mathrm{172}$) (Ito and Follows, 2005; Duteil et al., 2012):

$$\begin{array}{}\text{(4)}& \mathrm{\Delta }{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{Reg}}=r_{{\mathrm{P}:\mathrm{O}}_{\mathrm{2}}}\times \mathrm{\Delta }\mathrm{AOU}.\end{array}$$

AOU estimates the oxygen consumed during respiration and can be calculated as the difference of dissolved oxygen concentration (O₂) from the saturated concentration of oxygen (O₂^sat) (Weiss, 1970):

$$\begin{array}{}\text{(5)}& \mathrm{AOU}={{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{sat}}(T,S)-{\mathrm{O}}_{\mathrm{2}}.\end{array}$$

This can then be used to infer the efficiency of the biological pump (BP_Eff) as per Ito and Follows (2005):

$$\begin{array}{}\text{(6)}& \stackrel{\mathrm{‾}}{{\mathrm{BP}}_{\mathrm{Eff}}}=\stackrel{\mathrm{‾}}{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{Reg}}}/\stackrel{\mathrm{‾}}{{\mathrm{PO}}_{\mathrm{4}}},\end{array}$$

where $\stackrel{\mathrm{‾}}{X}$ is the mean value of X.

The contribution to DIC from the carbonate pump is estimated by

$$\begin{array}{}\text{(7)}& \mathrm{\Delta }{\mathrm{C}}_{{\mathrm{CaCO}}_{\mathrm{3}}}=\mathrm{0.5}\times \left[\mathrm{\Delta }\mathrm{ALK}+r_{\mathrm{N}:\mathrm{P}}\times \mathrm{\Delta }{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{Reg}}\right],\end{array}$$

where the term $r_{\mathrm{N}:\mathrm{P}}\times \mathrm{\Delta }$PO₄^Reg accounts for the reduction in ALK from production of nitrate (NO${}_{\mathrm{3}}^{-}$) (Sarmiento and Gruber, 2006), which is estimated using PO₄^Reg and the N:P stoichiometric ratio ($r_{\mathrm{N}:\mathrm{P}}=\mathrm{16}$). Finally, the preformed carbon concentration (C_pre) is obtained as follows:

$$\begin{array}{}\text{(8)}& \mathrm{\Delta }{\mathrm{C}}_{\mathrm{pre}}=\mathrm{\Delta }\mathrm{DIC}-\mathrm{\Delta }{\mathrm{C}}_{\mathrm{org}}-\mathrm{\Delta }{\mathrm{C}}_{{\mathrm{CaCO}}_{\mathrm{3}}}.\end{array}$$

3 Results

Changes in the climate system are presented first (Sect. 3.1), followed by their effects on air–sea gas exchange in the SO (Sect. 3.2). To understand these, we next quantify the different contributors to changes in surface pCO₂ (Sect. 3.3) and finish by analysing deep-ocean changes and the global oceanic carbon inventory (Sect. 3.4). Unless stated otherwise, the SO is defined as the ocean area south of 40^∘ S.

3.1 Ocean dynamics and sea ice cover

As a result of the insolation anomalies and associated feedbacks, the global mean annual SST anomaly at the LIG compared to PI as simulated by the ACCESS-ESM1.5 is equal to 0.17 ^∘C, with a pronounced warming at high latitudes and maximum positive SST anomalies over the North Atlantic (up to 4 ^∘C, Yeung et al., 2021). Our simulated temperatures are in line with the range of PMIP4 lig127k simulations (Otto-Bliesner et al., 2021). A model–data comparison of the LIG climate state is presented in Yeung et al. (2021).

A mean 0.53 ^∘C warming is simulated over the SO south of 40^∘ S (Fig. 1a). Warmer conditions are simulated everywhere south of 50^∘ S apart from a ∼ 1 ^∘C cooling centred at 58^∘ S in the South Atlantic and up to a 1.5 ^∘C cooling in the subantarctic eastern Pacific. The strongest warming is simulated over the Southeast Atlantic and Indian Ocean sectors (up to 3 ^∘C). Regional SSTs up to 4 ^∘C higher are simulated around 60^∘ S for both austral spring (Fig. A1b) and summer (not shown). The higher SSTs over the SO compared to PI are accompanied by a marked reduction in sea ice extent over both austral summer and winter (Fig. 1a), peaking at 41 % reduction in austral winter (Fig. A1a).

Figure 1 Annual mean (a) SST anomalies (^∘C) between LIG and PI overlaid with 15 % sea ice concentration (black for PI, magenta for LIG, solid lines for DJF, dashed lines for JJA). The value in the box shows the SO (40–90^∘ S) mean ΔSST. (b) Anomalies of the Ekman pumping velocities between LIG and PI (ms⁻¹) overlaid with sea ice concentration (as in a). (c) Global mean meridional streamfunction (Sv) for PI and (d) LIG. Positive values indicate clockwise water mass transport, and negative values indicate anticlockwise transport. Zonal mean (e) winds (ms⁻¹), and (f) Ekman pumping velocities (ms⁻¹) over the SO for PI (blue) and LIG (red).

A 1.5^∘ equatorward shift of the SH westerlies is simulated at the LIG (Fig. 1e), with a 10 % weakening of the winds south of 50^∘ S. This leads to ∼ 10 % weaker upwelling south of 55^∘ S and up to ∼ 20 % stronger upwelling north of 55^∘ S (Fig. 1f). Seasonally, the largest changes in upwelling are found for the winter and spring seasons (Fig. A1c and d). The equatorward shift of the westerlies reduces the northward Ekman transport south of 55^∘ S, leading to warming, while the higher Ekman transport north of 55^∘ S induces a cooling around 50^∘ S (Fig. 1a and b). The Antarctic Circumpolar Current is weaker, and the transport through Drake Passage is reduced from PI by ∼ 55 Sv (not shown).

3.2 Response of the air–sea gas exchange

These physical changes in the SO impact the carbon cycle and the air–sea CO₂ exchange (Fig. 2). It is worth reiterating here that we analyse changes in the carbon cycle using tracers for a constant atmospheric CO₂ concentration of 280 ppm for both simulations (Sect. 2.1), which enables us to solely analyse the impact of climatic and oceanic circulation changes on the carbon cycle. The CO₂ flux over the SO at PI shows a carbon uptake near the Antarctic coast, an outgassing band between 65–45^∘ S, and another uptake zone further north (Fig. 2a and blue line in d). This CO₂ outgassing is due to the upwelling of DIC-rich deep waters (Fig. 1f). At the LIG, the upwelling region widens over the sectors of the Atlantic and Indian oceans, and narrows over the eastern Pacific Ocean sector (Fig. 2). There is a strong increase in outgassing in the sectors of the Atlantic and western Indian oceans, as well as at ∼ 67^∘ S in the eastern Pacific Ocean sectors (Fig. 2c and red line in d). An increase in CO₂ uptake is simulated in the Amundsen, Bellinghausen, Weddell, Lazarev, Riiser-Larsen, and Ross seas and the subantarctic eastern Pacific (Fig. 2c). The zonal mean CO₂ flux in Fig. 2d shows a small increase in uptake south of 62^∘ S, possibly due to reduced mixing and reduced winter sea ice cover. Overall, there is a net increase in the mean SO outgassing by ∼ 150 % at the LIG compared to PI (Fig. 2e). This increase in outgassing mostly occurs during the austral winter and spring seasons (Fig. A1e and f).

Figure 2 Annual mean air–sea CO₂ flux (mol C m⁻² yr⁻¹) for (a) PI, (b) LIG, and (c) LIG−PI. Red colours indicate outgassing of CO₂ from the ocean, and blue colours indicate uptake by the ocean. Thick black lines show the zero line contour. (d) Zonal mean CO₂ flux (mol C m⁻² yr⁻¹) over the SO for PI (blue) and LIG (red). (e) Mean CO₂ flux over the SO (mol C m⁻² yr⁻¹) for PI (blue) and LIG (red).

This increased outgassing over the SO is compensated by increased uptake over other ocean basins, especially the North Pacific subpolar gyre, the South Atlantic Ocean, and the Northwest Indian Ocean, while higher outgassing is simulated in the North Atlantic (Fig. A2). In this paper, we focus primarily on the SO. To better understand the SO changes, we decompose the oceanic pCO₂ changes into their different components.

3.3 Changes in surface carbon dioxide partial pressure

The air–sea gas exchange is primarily controlled by the seawater partial pressure of CO₂ (pCO₂). Using the equations presented in Sect. 2.2, Fig. 3 shows a decomposition of the changes in pCO₂ into their SST, SSS, DIC, and ALK contributions. In line with the net increase in SO outgassing at the LIG (Fig. 2), the net pCO₂ over the SO is 1.2 µatm higher at the LIG compared to PI (black circle) (Fig. 3a). In agreement with the changes in air–sea CO₂ flux, surface pCO₂ is higher in a zonal band centred at ∼ 55^∘ S, while it is lower in coastal regions (contours in Fig. 3b and black line in Fig. 3g). The contributions from individual components add up to reflect the simulated differences reasonably well in terms of both magnitude (1.4 µatm, grey square in Fig. 3a) and spatial distribution (Fig. 3b), affirming the validity of the decomposition method. The overall slightly positive pCO₂ anomaly at the LIG compared to PI results from the competing effects of lower CO₂ solubility (red triangle, +5.65 µatm) and changes in surface DIC and ALK (blue diamond, −4.25 µatm).

Figure 3 Attribution of changes in annual mean surface pCO₂ over the SO (µatm). (a) Summary of decomposition of pCO₂ over the SO. Simulated difference in pCO₂ (black circle) and sum of contributions from all components (grey square, Sect. 2.2). pCO₂ change from solubility (SST + SSS, red triangle) and the sum of DIC and ALK components (blue diamond) are given in the same panel. Individual contributions from SST (brown triangle), SSS (magenta triangle), DIC (cyan diamond), and ALK (green diamond) are also shown. (b) Map of the sum of pCO₂ contributions from all components (corresponding to grey square in a) in shading, overlaid with the pCO₂ change simulated by the model as contours (corresponding to black circle in a). Maps of individual pCO₂ contributions from (c) SST (corresponding to the brown triangle in a), (d) DIC and ALK (corresponding to the blue diamond in a), (e) DIC (corresponding to the cyan diamond in a), and (f) ALK (corresponding to the green diamond in a). (g) Zonal mean contributions to pCO₂ over the SO (µatm) for simulated ΔpCO₂ (black), sum of contributions from all components (grey), solubility (red), sum of DIC and ALK components (blue), SST (dashed brown), SSS (dashed magenta), DIC (dashed cyan), and ALK changes (dashed green) are each shown individually. Note the non-linear scale given for the top and bottom thirds of (g).

The largest contributor to solubility is SST (brown triangle), while the largest contributor of the combined DIC and ALK effect is DIC (cyan diamond in Fig. 3a). Higher SSTs in the SO at the LIG lead to a 5.8 µatm pCO₂ increase (brown triangle in Fig. 3a). This SST-induced increase is present over most of the SO and is highest over the South Atlantic, Indian, and West Pacific regions (Fig. 3c). Changes in SSS do not contribute significantly to the pCO₂ anomalies (magenta triangle in Figs. 3a and A3d). Changes in DIC cause the largest single contribution to the overall pCO₂ change, with the net decrease in surface DIC leading to a −10.9 µatm pCO₂ change (cyan diamond in Fig. 3a). As an exception to this, the higher DIC concentrations in the Atlantic Ocean around 55^∘ S and in the eastern Pacific sector between 60–70^∘ S (Figs. A1g, h, and A4d) lead to higher pCO₂ (Fig. 3e). This higher DIC is due to increased upwelling in these regions (Figs. 1b and A1c, d), resulting from the equatorward shift of the SH westerlies. The contributions based on changes in ALK and DIC have very similar patterns, albeit of opposite signs, with a small difference over the West Pacific around 50^∘ S. The combined effect of ALK and DIC changes leads to a total net decrease in pCO₂ by 4.25 µatm due to the higher impact of changes in DIC (−10.9 µatm, cyan diamond) compared to ALK (+6.65 µatm, green diamond). Overall, south of 45^∘ S, changes in SST lead to an overall increase in pCO₂, with the highest contributions over the Indian and West Pacific sectors of the SO, while changes in DIC lead to an overall decrease of pCO₂, with the strongest reductions in the South Atlantic, Indian, and West Pacific sectors of the SO but an increase in the Amundsen Sea sector (Fig. 3).

The zonal mean patterns of these contributions are presented in Fig. 3g. Coastal regions around Antarctica south of ∼ 75^∘ S show that the lower ΔpCO₂, both simulated (black line) and calculated (grey line), primarily arise from the combined DIC and ALK component (blue line) that is mainly attributed to the lower DIC (dashed cyan line) compensated by changes in ALK (dashed green line). Between 70 and 60^∘ S, ΔpCO₂ anomalies are close to zero as the reduced CO₂ solubility due to higher SSTs (red line) is compensated by the combined DIC and ALK changes (blue line). Between 60 and 45^∘ S, the simulated surface ocean ΔpCO₂ is positive (∼ 2 µatm, black line) as the effect of higher SSTs dominates over the combined ALK and DIC components. North of ∼ 45^∘ S, the positive ΔpCO₂ signal arises from the DIC + ALK component, while solubility mitigates the anomaly (Fig. 3g).

Figure 3g shows that while solubility is mostly controlled by changes in SST, contributions from SSS changes are important south of 72^∘ S. Similarly, DIC has the dominant control over the DIC + ALK component. The contributions from both the solubility and (DIC + ALK) components reach their maxima (solubility being positive and DIC + ALK being negative) around 62^∘ S, which corresponds to the maximum divergence in wind-driven surface currents (Fig. 1f). To summarize, the higher overall pCO₂ (and hence lower air–sea CO₂ flux) at the LIG compared to PI can be attributed to higher SSTs south of 45^∘ S, and the combined DIC + ALK component between 45 and 35^∘ S due to the equatorward shift of the upwelling regions (Fig. 3g). The SST patterns have already been discussed in Sect. 3.1. DIC patterns (Fig. A4d) can result from changes in circulation and the biological pump. These are investigated in the next section.

Figure 4 Global marine carbon budget and decomposition. Global zonal and annual mean anomalies of (a) dissolved oxygen concentration (µmol kg⁻¹), (c) phosphate concentration (µmol kg⁻¹), (e) regenerated phosphate concentration (µmol kg⁻¹), (g) preformed phosphate concentration (µmol kg⁻¹), (b) total DIC (µmol kg⁻¹), (d) remineralized carbon (µmol kg⁻¹), (f) dissolved carbonate (µmol kg⁻¹), and (h) preformed DIC (µmol kg⁻¹) are given in the individual panels. Note that the phosphate components (c, e, g) and the DIC components (b, d, f, h) each share common colour bars.

3.4 Deep-ocean changes and carbon inventory

Figure 1c and d show that North Atlantic Deep Water (NADW) formation is ∼5 Sv higher in our LIG simulation compared to the PI simulation, leading to colder waters at depths of 1–2 km between 0–60^∘ N, accompanied by higher SSTs and SSSs in the Labrador Sea (not shown). We also simulate a ∼ 4 Sv reduction in AABW formation (Fig. 1c and d). This reduction in AABW formation results in a warming of SO deep waters by up to 2 ^∘C compared to PI. Due to a northward shift and weakening of winds (Fig. 1e, f), the Antarctic Intermediate Water (AAIW) formation regions shift northward and the AAIW formation rate slows down (Downes et al., 2017). The higher SST, equatorward-shifted upwelling, and reduced sea ice extent lead to an increase in net primary production (NPP) and export production over the SO of ∼ 17 % and ∼ 11 % respectively (Fig. A4a and b).

These circulation changes impact the DIC distribution in the ocean (Fig. 4b). For instance, the reduced formation rate of AAIW (Fig. 1d) increases residence times and leads to lower dissolved oxygen (Fig. 4a), higher PO₄ (Fig. 4c), higher remineralized carbon (Fig. 4d), and higher DIC concentrations (Fig. 4b) at intermediate depths of the SO north of 55^∘ S. Similarly, the increased accumulation of nutrients and DIC, as well as oxygen depletion in deep and abyssal waters (Fig. 4a, b and c), can be attributed to a weaker AABW formation rate (Fig. 1d), which results in a higher efficiency of the biological pump, as detailed below.

The ACCESS-ESM1.5 simulates a ≥50 µmol kg⁻¹ decrease in oxygen everywhere below 3 km at the LIG compared to PI. A decrease in oxygen can also be seen in the SO across all depths, including a northward-reaching tongue attributable to AAIW at ∼ 1 km depth extending to the Equator (Fig. 4a). DIC anomalies follow the patterns of dissolved oxygen, albeit with the opposite sign. The simulated DIC concentrations are ∼ 50 µmol kg⁻¹ higher across all depths in the SO and globally below 3.5 km depth (Fig. 4b).

DIC anomalies are decomposed into contributions from remineralized organic carbon, dissolved calcium carbonate, and preformed carbon components (Fig. 4d, f and h, Sect. 2.3). A total of 60 % of the DIC increase in the SO and abyssal ocean can be attributed to an increase in remineralized carbon (Fig. 4d) resulting from a 10 % more efficient biological pump (Eq. 6, Sect. 2.3). This increase in remineralized carbon can be attributed to increased residence times, which are in turn due to weaker bottom- and intermediate-water formation rates (Fig. 1c and d). Weaker AABW and AAIW indeed lead to positive apparent oxygen utilization (AOU) and regenerated PO₄ anomalies (Fig. 4e). In addition, a 2 ^∘C warming of bottom waters leads to a ∼ 6 % increase in remineralization rates, thus further contributing to the higher remineralized carbon. Around 25 % of the abyssal increase in DIC is attributed to the carbonate pump (Fig. 4f). This reflects changes in NPP, export production (Fig. A4), and water mass residence times, given that carbonate production is a constant percentage of total NPP in this model setup and carbonate dissolution is constant and independent of water chemistry or temperature. Preformed DIC represents only ∼ 5 % of the DIC changes (Fig. 4h). This increased sequestration of DIC in the deep ocean reduces the surface DIC concentration (Figs. A4 and A1g, h), thus contributing to a lowering of surface pCO₂ (Fig. 3e and g).

For the Northern Hemisphere, stronger NADW (Fig. 1c and d) results in decreased DIC concentrations in intermediate and deep waters of the North Atlantic by up to 75 µmol kg⁻¹ (Fig. 4b). This negative DIC anomaly can be explained by a decrease in both remineralized and preformed DIC. Stronger NADW subducts more DIC-depleted surface waters into the deep ocean (Duteil et al., 2012) and reduces residence times, thus leading to a decrease in remineralized carbon. This is also associated with a ∼ 40 µmol kg⁻¹ increase in dissolved oxygen content (Fig. 4a) and an up to 0.75 µmol kg⁻¹ decrease in PO₄ concentration (Fig. 4c). The reduction in preformed carbon (Fig. 4h) is most likely due to a new NADW formation site in the Labrador Sea and a slight northward shift of the deep-water formation regions in the Norwegian and Greenland seas.

4 Discussion and conclusion

We analyse the SO marine carbon cycle response to warmer conditions as simulated in an equilibrium simulation of the LIG. The lig127k simulation performed with the ACCESS-ESM1.5 presented here displays an annual mean warming of 0.53 ^∘C at the surface of the SO compared to PI. This simulated southern high-latitude warming (Fig. 1a) is in agreement with the multi-model mean of the PMIP4 lig127k simulations, although it is at the higher end of the spectrum (Otto-Bliesner et al., 2021). Seasonally, we simulate regional SSTs up to 4 ^∘C higher around 60^∘ S for both the austral spring and summer, in line with both proxy records (Capron et al., 2017; Hoffman et al., 2017) and the PMIP3 125 ka experiment performed with the NORESM-1ME model (Kessler et al., 2018), while the PMIP4 lig127k multi-model mean for austral summer displays less than 0.5 ^∘C warming over the SO (Otto-Bliesner et al., 2021). The simulated SO sea ice concentration at the LIG for austral winter (Figs. 1 and A1) is also in good agreement with those inferred in previous studies (Holloway et al., 2016, 2017), even though the changes in Antarctic sea ice area of the ACCESS-ESM1.5 are the largest amongst all the CMIP6-PMIP4 models (Otto-Bliesner et al., 2021).

A 1.5^∘ equatorward shift of the westerlies, resulting in a 10 % weakening of the westerlies south of 50^∘ S, is simulated at the LIG (Fig. 1e), leading to significant changes in SO upwelling (Fig. 1b and f). There is no clear consensus for position and strength of the SH westerlies during the LIG (Fogwill et al., 2014), although the higher obliquity might have led to a weakening of the westerlies compared to PI (Timmermann et al., 2014). These changes in winds lead to a northward shift of the AAIW formation regions. Our LIG simulation is also characterized by a weakening of AABW formation rates, which might be due to changes in surface density. Significant weakening of AABW during the LIG due to Antarctic meltwater discharge has previously been inferred from paleo-proxy records (Hayes et al., 2014; Rohling et al., 2019). A shift in westerlies might also contribute to a weakening of the AABW (Menviel et al., 2008; Huiskamp et al., 2016; Glasscock et al., 2020).

The reduced upwelling south of 55^∘ S and weaker AABW transport simulated here lead to an increased sequestration of DIC in the deep ocean, through an increase in the efficiency of the biological pump (Fig. 4). This reduces the surface DIC concentrations, leading to a net reduction in outgassing over the SO, as has been previously hypothesized (Toggweiler et al., 2006). However, while the combined DIC and ALK contributions would lead to a lower pCO₂ at the surface of the SO, this change is overcompensated by the warmer conditions (Fig. 3). Reduced solubility due to higher SSTs leads to an increase in outgassing over most of the SO, while the reduced sea ice cover does not seem to significantly impact the CO₂ fluxes. Assessing the impact of sea ice changes on CO₂ fluxes (Appendix A), we find that reduced sea ice concentration at the LIG in the Weddell and Ross seas leads to a 5 % increase in CO₂ uptake in autumn and winter (Fig. A5). Reduced solubility results in a net outgassing of CO₂ over the SO, with a ∼150% increase at the LIG compared to PI (Fig. 2), and the largest increase over the austral winter and spring seasons (Fig. A1). The simulated NPP and export production are ∼17% and ∼ 11 % higher, respectively, over the SO at the LIG compared to PI (Fig. A4), providing a negative feedback to the higher outgassing. Although the simulated nutrient and alkalinity concentrations are underestimated over the SO, the vertical gradients of these are captured reasonably well (Ziehn et al., 2020). As such, these biases should not significantly impact our results. All of our analysis is based on a constant atmospheric CO₂ concentration of 280 ppm to allow quantification of the effects of the LIG climate on the carbon cycle independently of the background CO₂ concentration. However, this constant atmospheric CO₂ concentration neglects feedbacks related to CO₂ uptake and outgassing. Nevertheless, the lower CO₂ at LIG (275 ppm) compared to PI (284.3 ppm) would suggest the LIG SO to be an even greater CO₂ source to the atmosphere, implying a stronger sink somewhere else in the ocean or on land (Brovkin et al., 2016).

Numerical studies have unequivocally shown the impact of changes in the magnitude of the westerlies on the carbon cycle, with stronger westerlies leading to increased upwelling and SO outgassing and vice-versa (e.g. Menviel et al., 2008, 2018; d'Orgeville et al., 2010; Lauderdale et al., 2013, 2017; Huiskamp et al., 2016; Gottschalk et al., 2019); however, the impact of changes in the latitudinal position of the westerlies on the carbon cycle is less certain (e.g. d'Orgeville et al., 2010; Völker and Köhler, 2013; Lauderdale et al., 2013; Huiskamp et al., 2016; Gottschalk et al., 2019). Given the simulated increase in DIC content in the deeper SO and the reduced DIC at the surface of the SO, our simulations suggest that the simulated equatorward shift of the westerlies reduces the SO CO₂ outgassing.

SH westerlies are projected to strengthen and shift poleward over the coming century (Collins et al., 2013; Zheng et al., 2013; Goyal et al., 2021), contrary to the LIG simulations presented here. An increase in SO NPP at the LIG is simulated here, in line with CMIP6 projections of the coming century (Kwiatkowski et al., 2020). Thus, changes in the carbon cycle simulated at the LIG may not serve as a good analogue for potential future changes. Nevertheless, the simulated enhanced SO CO₂ outgassing, despite a slight equatorward shift of the westerlies, supports a weaker capability of the Southern Ocean to take up anthropogenic CO₂ over the coming century.

Appendix A: Effect of sea ice on air–sea gas exchange

To estimate the effect of sea ice changes at the LIG on the air–sea gas exchange, we use a modified version of the equation of CO₂ flux from Wanninkhof (2014):

$$\begin{array}{}\text{(A1)}& F_{{\mathrm{CO}}_{\mathrm{2}}}=\mathrm{7.7}\times {\mathrm{10}}^{-\mathrm{4}}\times \left|U^{\mathrm{2}}\right|\times \mathrm{\Delta }{\mathrm{pCO}}_{\mathrm{2}}\times (\mathrm{1}-\mathrm{Sc}),\end{array}$$

where $F_{{\mathrm{CO}}_{\mathrm{2}}}$ is the CO₂ flux (mol C m⁻² yr⁻¹), U is the average wind speed (m s⁻¹), ΔpCO₂ is the difference in partial pressure of CO₂ between ocean surface and atmosphere (µatm), and Sc is the sea ice concentration. To investigate the effect of LIG sea ice changes on the LIG CO₂ flux, we estimate $F_{{\mathrm{CO}}_{\mathrm{2}}}$ using both LIG Sc and PI Sc, while keeping all the other variables at LIG values. The resulting patterns are presented in Fig. A5. Figure A5 shows that the reduced sea ice during LIG compared to PI leads to a less than 5 % increase in carbon sink over the Ross Sea region all year round and over the Weddell Sea region in autumn and winter and an around 2 % increase in outgassing in the Lazarev and Weddell seas.

Figure A1 Seasonal anomalies for austral winter, (JJA, a, c, e, g) and austral spring (SON, b, d, f, h) of (a, b) SST (^∘C) overlaid with 50 % sea ice concentration (dashed for PI and solid for LIG), (c, d) Ekman pumping velocities (ms⁻¹) overlaid with 50 % sea ice concentration (dashed for PI and solid for LIG), (e, f) sea–air CO₂ flux (mol C m⁻² yr⁻¹, with positive indicating outgassing from and negative indicating uptake by the ocean), and (g, h) surface DIC concentrations (µmol kg⁻¹).

Figure A2 (a) Annual mean air–sea CO₂ flux (mol C m⁻² yr⁻¹) for LIG−PI. Red colours indicate outgassing of CO₂ from the ocean, and blue colours indicate uptake by the ocean. Thick black lines show the zero line contour. (b) Zonal mean CO₂ flux (mol C m⁻² yr⁻¹) for PI (blue) and LIG (red).

Figure A3 Maps of individual pCO₂ contributions (µatm) from (a) SST (corresponding to the brown triangle in Fig. 3a) and (b) SSS (corresponding to the magenta triangle in Fig. 3a).

Figure A4 Changes in annual mean productivity and related variables. (a) Anomalies of depth-integrated annual mean NPP (as chlorophyll in seawater, kg m⁻³) overlaid with 15 % sea ice concentration (black for PI, red for LIG, solid lines for DJF, and dashed lines for JJA). (b) Anomalies of annual mean detrital concentration at 200 m depth (kg m⁻³) overlaid with sea ice (same as a). (c) Anomalies of annual mean surface phosphate concentrations (µmol kg⁻¹) overlaid with sea ice (same as a). (d) Anomalies of annual mean surface DIC concentrations (µmol kg⁻¹) overlaid with sea ice (same as a).

Figure A5 Calculated seasonal CO₂ fluxes using LIG ΔpCO₂, winds and (left) LIG sea-ice concentration, (middle) PI sea ice concentrations, and (right) differences in CO₂ flux for calculations using the LIG sea ice concentrations compared to PI sea ice concentrations. Note the reduced colour scale in the right-hand column. Calculations are based on Wanninkhof (1992) and Wanninkhof (2014). The red and blue contours indicate 15 % sea ice concentrations for LIG and PI, respectively.