Stratification constrains future heat and carbon uptake in the Southern Ocean between 30°S and 55°S

The Southern Ocean between 30°S and 55°S is a major sink of excess heat and anthropogenic carbon, but model projections of these sinks remain highly uncertain. Reducing such uncertainties is required to effectively guide the development of climate mitigation policies for meeting the ambitious climate targets of the Paris Agreement. Here, we show that the large spread in the projections of future excess heat uptake efficiency and cumulative anthropogenic carbon uptake in this region are strongly linked to the models’ contemporary stratification. This relationship is robust across two generations of Earth system models and is used to reduce the uncertainty of future estimates of the cumulative anthropogenic carbon uptake by up to 53% and the excess heat uptake efficiency by 28%. Our results highlight that, for this region, an improved representation of stratification in Earth system models is key to constrain future carbon budgets and climate change projections.

T he Southern Ocean is a dynamically complex region. The strong wind-driven Antarctic Circumpolar Current drives a residual overturning circulation, consisting of an upwelling of circumpolar deep waters around the Polar Front, a residual northward transport with gradual water mass transformation to Antarctic Intermediate and Mode Waters (IW and MW), and finally subduction under subtropical waters. The upwelled water mass is cold and undersaturated with respect to anthropogenic carbon, allowing it to efficiently absorb large amounts of atmospheric excess heat and anthropogenic carbon [1][2][3] . Therefore, the subduction of IW and MW masses, occurring approximately between 30°S and 55°S, provides one of the major gateways carrying anthropogenic carbon (C ant ) and excess heat (H excess ) into the interior ocean (e.g. refs. [4][5][6][7], where they stay isolated from the atmosphere on decadal to millennial timescales 8,9 . For the historical period from 1850 to 2005, it has been estimated that 43% of C ant and 75% of H excess have entered the ocean south of 30°S 10 , although the Southern Ocean accounts for only 30% of the total ocean surface area. Over the same period, the region between 30°S and 55°S is responsible for 27 and 50% of the global ocean C ant and H excess uptake despite covering only 21% of the world ocean, according to the models analysed in this study. Future projections of C ant and H excess uptake in the region between 30°S and 55°S from the last two generations Earth system models (ESMs) remain uncertain 11,12 (Fig. 1) because ESMs struggle to capture the complex dynamical and biogeochemical processes in this region [13][14][15] . Despite improvements in model performance in successive phases of the Coupled Model Intercomparison Project (CMIP), this progress might be too slow to warrant significantly reduced uncertainty of ESM projections within the next decade 16 . Since this is the time horizon for framing climate mitigation policies that allow for meeting stringent climate targets 17,18 , more efforts have to be put into model analysis, i.e., understanding the roots of this uncertainty and reducing uncertainty in key climate metrics such as the projected carbon and heat uptake. The technique of emergent constraints provides a means to constrain a model ensemble through an emergent strong statistical relationship between an observable quantity of current climate and future changes in a variable of interest 19,20 . It has been used to constrain several aspects of the terrestrial [21][22][23] and marine [24][25][26] carbon cycle.
In this work, we identify a key mechanism that explains the large inter model-uncertainty in future projections of H excess and C ant uptake between 30°S and 55°S across both CMIP5 and CMIP6 ESMs. We focus on the region between 30°S and 55°S since this is the area where intermediate and mode waters are formed and subducted (Methods). We find that the climatological stratification state in this region is tightly related to this subduction and we use this finding to robustly constrain both future H excess and C ant uptake. We note that our definition of C ant and H excess includes changes induced by climate change such as changes in ocean circulation, wind conditions and primary production (Methods).

Results
Linking stratification to oceanic C ant and H excess uptake. We find that stratification biases in CMIP5 and CMIP6 ESMs in the region between 30°S and 55°S are strongly related to the amount of their future uptake of excess heat per degree of transient global warming (H excess uptake efficiency) and anthropogenic carbon (Fig. 2). Models showing a positive density bias that increases with depth relative to the surface bias (indicating stronger-thanobserved stratification) tend to simulate a low uptake of C ant and low H excess uptake efficiency. The opposite is true for models that show an increasingly negative density bias profile relative to their surface bias (indicating weaker-than-observed stratification).
In order to develop an emergent constraint from this apparent relationship, we need to capture the characteristics of the vertical structure of these density profiles in one metric. To this end, we use a stratification index 27 , which is the cumulative sum of density differences with respect to surface density (Methods), here applied over the upper 2000 m of the water column. This depth range has been chosen as it encompasses the MW and IW formation and subduction pathways in CMIP ESMs 7,28 , and modern observational coverage is good, since it is covered by standard ARGO floats.
We identify the core of IW in each ESM by determining the depth of the salinity minimum at 30°S (ref. 28 ), and we find that the stratification index is highly correlated to both (1) the depth at which IW are subducted (R = −0.83) and (2) the subducted volume of IW and overlying MW (R = −0.77, both shown in Supplementary Fig. 1). Therefore, consistent with a previous study 29 , we find that the modelled volume of MW and IW formation is of high importance for determining the efficiency of C ant and H excess sequestration. However, the stratification index has the clear advantage of being straightforward to estimate from model output while the identification of water masses is more challenging and model-dependent 28 . We find that ESMs with high-stratification index and correspondingly low C ant uptake typically simulate lower uptake in the region around 55°S, but more importantly, the northward extent of their uptake is much more limited compared to models with low stratification index (Fig. 3). The latter models project accumulated uptake of more than 100 mol C m −2 in large regions north of 40°S in Pacific, Indian and Atlantic sectors, where low C ant -uptake models show uptake below 50 mol C m −2 (see also Supplementary Fig. 2 for a zonal mean view of C ant uptake and Supplementary Fig. 3 for an equivalent to Fig. 3 but for CMIP6 models). The higher (lower) C ant uptake simulated by ESMs with low (high) stratification index is connected to a steeper (shallower) surface-to-depth gradient of the anthropogenic component of dissolved inorganic carbon (DIC ant ) concentration along the vertical zonal mean section between 30°S and 55°S ( Fig. 3c, d). A similar surface-to-depth feature can also be seen in the warming efficiency (Fig. 3e, f).
It is physically plausible for a model that exhibits a stronger stratification than observed to take up less carbon and heat than a model with weaker than observed stratification. However, the fact that present-day stratification is related to projected future uptakes across our model ensemble is not obvious, since stratification is changing with progressing climate change. We find that, in the region between 30°S and 55°S, the projected stratification bias of models relative to each other remains largely unchanged, i.e., a model that simulates a stronger contemporary stratification than the multi-model mean will do so for future time periods, too. The correlation between the mean present-day (1986-2005) stratification index and the mean future (2080-2099) stratification index is 0.91 ( Supplementary Fig. 4).
The target variables for our constraints are (i) the cumulative ocean C ant uptake [Pg C] as this minimises interannual and decadal variability (compared to annual carbon fluxes) while preserving trends and (ii) the 20-year average of ocean H excess uptake efficiency [TW°C −1 ] defined as the ratio between ocean H excess uptake rate [TW] and global atmospheric surface warming [°C] (ref. 30 ). The latter choice is motivated by the fact that the simulated atmospheric surface temperature that forces the oceanic heat uptake rates depends on each model's response to radiative forcing. We, therefore, normalise the excess heat uptake by global surface warming [31][32][33] . Additional information on the validity of our findings for the H excess uptake rate [TW] without normalisation are presented in the supplement ( Supplementary  Fig. 5).
In our analysis of H excess uptake efficiency, we merge the CMIP5 and CMIP6 ensembles because both rely on scenarios with the same end of century radiative forcing (RCP8.5 and SSP5-8.5, respectively). Such a merger is not meaningful for the C ant uptake as the CMIP6 SSP5-8.5 scenario reaches considerably higher end-of-century atmospheric CO 2 concentrations than the CMIP5 RCP8.5 scenario 34 . The radiative forcing due to higher CO 2 concentrations in SSP5-8.5 is compensated by lower concentrations of other greenhouse gases, mainly methane and nitrous oxide.
Reducing uncertainties in ocean uptake projections. Significant negative correlations exist between the simulated present-day water-column stratification index and both cumulative C ant uptake and H excess uptake efficiency at the end of the century (Fig. 4a-c). We note that the correlations between the stratification index and H excess are still significant but less robust without normalising the H excess uptake ( Supplementary Fig. 5). These correlations indicate that a more stratified ocean absorbs less C ant and H excess . For C ant uptake, the high correlation with contemporary stratification is very stable over time (Fig. 4g), whereas for H excess uptake efficiency the correlation is initially low but gets stronger with time (see Discussion) and reaches values of 0.7 (P = 0.003) for CMIP5 and 0.8 (P < 0.001) for CMIP6 at the end of the century (Fig. 4h). The WOA13-based stratification index and its uncertainty are estimated as 64.08 ± 0.58 kg m −3 (Methods). This value is close to the CMIP5 and CMIP6 ensemble mean of 64.72 ± 3.80 kg m −3 and 65.02 ± 1.70 kg m −3 , respectively. However, the model spread around the mean is substantial for both CMIP5 and CMIP6, with CMIP5 having a model uncertainty that is more than twice as large as the one of CMIP6.
Based on the high correlations for both the CMIP5 and CMIP6 ensembles, we apply the emergent constraint approach (Methods) to constrain the uncertainty in future projections of C ant uptake and H excess uptake efficiency. We assume that all models are independent as done in other studies 35 . We note that this is a limitation of our study as some of these ESMs share components and code. Likewise, many CMIP6 models have been developed starting from their predecessor CMIP5 models such that the two ensembles are not entirely independent. An alternative approach could be based on an adaptive model weighting scheme or an ensemble reduction 36,37 , but this is beyond the scope of our study. After applying the observational constraint (Methods), the uncertainties of the cumulative C ant uptake between 30°S and 55°S are considerably reduced by 53 and 32% for CMIP5 and CMIP6, respectively. The associated best estimate of cumulative C ant uptake increases by 3 and 6% for CMIP5 and CMIP6 respectively compared to the prior-constraint estimate (Table 1). Similarly, the after-constraint uncertainty of H excess uptake efficiency for the combined CMIP5/CMIP6 ensemble is strongly reduced by 28% and the associated estimate increases by 7%.

Discussion
Our emergent constraint identifies a strong link between contemporary stratification in CMIP5/6 models and their ability to continuously take up C ant and H excess under a high-CO 2 future scenario in the Southern Ocean between 30°S and 55°S. The ESMs' stratification index correlates strongly with (i) the simulated depth at which the IW (and the overlying MW) are subducted and (ii) the simulated subducted water volume, here loosely referred to as the volume above the IW core (both shown in Supplementary Fig. 1). This suggests that a deeper position of the IW core is accompanied by a larger subduction volume, and hence a more efficient C ant and H excess sequestration in our model ensemble. This importance of the volume of ventilated waters for future C ant uptake in the Southern Ocean has been found in an independent study 29 . We note that the relationship between contemporary stratification and C ant and H excess uptake worsens . c-f Corresponding zonally-averaged transects from the c, e low-and d, f high-stratification models for dissolved inorganic C ant concentration and warming efficiency. The latter is defined as the ratio between the transient minus piControl temperature anomaly and the global surface atmospheric warming. The black dotted lines in panels c-f represent the zonally averaged density isolines crossing the depth of the salinity minimum at 30°S 28 . See Supplementary Fig. 2 for the corresponding CMIP6 models.  Fig. 4 Emergent constraints on the sensitivity of projected carbon and heat uptake to stratification. a, b cumulative future C ant uptake between 30°S and 55°S and c future ocean H excess uptake efficiency and their emergent relation to the contemporary stratification. CMIP5, CMIP6 and combined CMIP5/ 6 ensembles are denoted in blue, red and magenta, respectively. Emergent constraints include a linear regression fit (dotted line), its 68% prediction interval (abbreviated pred. int., coloured shaded area according to the ensemble), the observational constraint (vertical black dots) and its uncertainty (grey shaded area). All emergent constraints are accompanied by d-f prior-and after-constraint probability density functions and g, h correlation timeseries obtained by sliding the predictand over time while leaving the predictor fixed. when extending the region of interest south of 55°S and outside of the area of IW and MW subduction. Here, the C ant and H excess uptake is sensitive to other processes, such as sea-ice dynamics and bottom-water formation at the southward limb of Southern Ocean overturning circulation, and we find no direct link to the contemporary stratification ( Supplementary Fig. 6). It has been shown before that formation of mode and intermediate waters is key for carbon and heat uptake 5,38 , and that this water mass formation appears to be linked to the simulated winter mixed layer depths of CMIP5 models 14 . We find, however, that the relationship between future cumulative C ant and H excess uptake efficiency and mixed layer depth in our region is only weak. There is no significant correlation between annual mean or maximum winter mixed layer depth and carbon and heat uptakes ( Supplementary Figs. 7, 8). Deep winter mixing in our region is thought to be the main contributor to carbon and heat subduction, and therefore such relatively low correlations might seem surprising. A previous study 14 has shown that stratification biases in ESMs contribute to setting the maximum winter mixed layer depth, and this effect is also captured by our constraint. In addition, other processes (such as diapycnal mixing and Ekman pumping) are also important for the eventual subduction of carbon and heat away from the seasonally varying base of the mixed layer 39,40 . A robust relationship can be found when taking the stratification of the upper 2000 m of the water column into account, but we note that our constraint is not sensitive to the exact lower bound of this depth range when it is varied between 1000 and 2000 m. We translate these findings into a physically plausible and robust emergent constraint for future projections of two generations of ESMs. It provides us with the unique possibility to constrain the model uncertainty for two highly important quantities at the same time.
For the CMIP5/6 generation of models, we find that the simulated contemporary stratification between 30°S and 55°S is highly correlated to its projected future values across CMIP5/6 (R correlation of 0.91). Models with a strong stratification store more of the excess heat in the upper ocean, thereby creating stronger stratification changes, while the opposite is true for weakly stratified models 41 . This mechanistic explanation suggests that ESMs that simulate a realistic contemporary density profile are more reliable in simulating future density profiles. The predictor of our emergent constraint, i.e. the contemporary stratification, hence also constrains the future stratification and, as demonstrated here, the future cumulative C ant uptake and the H excess uptake efficiency. A recent study indicates that projected patterns of heat storage are primarily dictated by the preindustrial ocean circulation 42 . Contemporary oceanic storage of anthropogenic carbon and excess heat have distinct patterns [43][44][45] and redistributed heat and carbon are projected to have opposing signs, leading to a more horizontal structure of heat storage than seen in the patterns of carbon storage in the Southern Ocean 46 . However, these differences are reduced in future projections of the Southern Ocean as the spatially averaged added heat becomes dominant over the redistribution of heat due to circulation changes 42,46 . This is consistent with our findings, specifically with the initially low but increasing correlation between stratification and excess heat uptake efficiency (Fig. 4e). We note that other quantities like the nutrient cycle and primary productivity are also closely linked to stratification, e.g. it has been shown that CMIP5 models with a stronger bias in contemporary surface stratification tend to predict larger climate-induced declines in surface nutrients and net primary production 47 . Due to the high importance of contemporary stratification biases for future marine projections, it is essential to reduce them.
Our results identify significant stratification biases for most CMIP5/6 models in the areas of MW and IW formation, but also that the representation of stratification in the Southern Ocean between 30°S and 55°S has improved between CMIP5 and CMIP6. In fully coupled ESMs, it remains difficult to identify the ultimate source of biases. The emergent-constraint method is only able to identify systematic biases associated with the variables used in the emergent-constraint relationships. It does not highlight missing processes or dynamical biases common in ESMs which are not directly related to the observable processes or variables used in the constrained process 35,37,48 . As many of the model biases in the Southern Ocean temperature and salinity structure are concentrated in recently ventilated layers or in the deep Atlantic, they appear to stem from inaccuracies in the North Atlantic Deep Water formation regions or in the surface climate over the Southern Ocean 16 . Recently, a strong emergent constraint relationship has been found between surface salinity and cumulative C ant uptake in the Southern Ocean 49 . In combination with our constraint, this indicates that surface salinity is a fundamental player setting the Southern Ocean stratification. Here, the upper ocean properties like salinity are highly sensitive to a multitude of uncertainties in the sea-ice, ocean, and atmosphere components of an ESM, e.g. westerly jet position, Antarctic seaice extent and its potential relation to precipitation, clouds, mixing and transport by eddies. It takes a tremendous effort to model or parameterise all these processes in a realistic manner and a significant reduction of bias is not to be expected within this decade 16 . For the ocean, it has been found that eddy-induced diffusion is an important factor in setting the simulated stratification 41 . Hence, a better representation of eddies, be it through increased eddy-resolving resolution 50 or through improved eddy parameterisations 51 will very likely contribute to reducing stratification biases in ESMs. Future studies should elucidate processes that could contribute to the bias and large spread in the stratification index simulated across ESMs, for instance, our stratification index could be influenced by the water mass properties of the circumpolar deep water, which is formed in the North Atlantic. A better understanding of the linkage between North Atlantic climate representation and the Southern-Ocean water-mass properties across ESMs could be valuable.
Ensembles of ESMs remain our only tool at hand to investigate the response of the Earth system to future scenarios of anthropogenic forcing. Reducing the large uncertainties arising from, among others, the representation of Southern Ocean dynamics in these models remains a challenge. The identification of emergent constraints, such as the one presented here, are invaluable as they can help to guide model development and, importantly, to speed up the provision of critical knowledge on expected future changes.  28 . We use a single ensemble member (r1i1p1(f1) or equivalent) per model. The selected ESMs provide full periods of the following three standard CMIP5 (CMIP6) experiments: piControl, historical and RCP8.5 (SSP5-8.5). For our study, ocean H excess uptake, ocean C ant uptake and global atmospheric surface warming are calculated using the air-sea heat flux, the air-sea CO 2 flux and the surface air temperature, respectively. The anthropogenic or excess component is obtained as the difference between the historical or future scenario and the preindustrial control experiments. Thus, C ant and H excess include changes induced by climate change (e.g. changes in ocean circulation, wind conditions, primary production).

Methods
Latitudinal extent of the region considered for the constraint. In our study, we focus on the area where intermediate and mode waters are formed and subducted as these processes are the main drivers of ocean carbon and heat uptake in the Southern Ocean 10,28,38 . This area lies between 30°S and 55°S in all CMIP5 and CMIP6 models. The 30°S is a commonly used northern boundary for the Southern Ocean and its subduction region 5,10,21 . The 55°S southern boundary is chosen to exclude the influence of sea-ice on air-sea fluxes, which would complicate the uptake-stratification relationship (Supplementary Fig. 9). According to the CMIP5 and CMIP6 zonal wind stress distribution 16 , the 55°S southern boundary generally excludes the southward limb of the Southern Ocean overturning circulation, that is not related to the subduction process of interest in this study.
Density calculations and stratification index. We calculated in situ density (⍴) from each ESM´s potential temperature and practical salinity (after conversion to absolute salinity and conservative temperature) following TEOS-10 standards 53 . Three-dimensional ⍴ fields have been area-weighted and averaged along horizontal surfaces to produce one-dimensional vertical profiles in native (model-dependent) vertical resolution.
We use a Stratification Index (SI) based on ref. 27 to characterise the stratification of the water column: where z 0 is the sea surface and z i = z i−1 + 200 for i = 1, …, 10 Probability density functions for the emergent constraints. The prior probability density functions for cumulative C ant uptake and H excess uptake efficiency assume that all models are equally likely to be correct and lead to a Gaussian distribution, and so is the probability density function of the observational constraint P(x) 21 . The probability density function of the constrained estimate P(y) was generated following established methodologies by normalising the product of the conditional probability density function of the emergent relationship P(y│x) and the probability density function of the observational constraint P(x): 21,22,25,26 where x and y are the predictor and the predictand, respectively. σ f = σ f (x) and is the 'prediction error' of the emergent linear regression.
Observational constraint. The World Ocean Atlas 2013 version 2 (WOA13) annual climatology of ⍴ (refs. 54,55 ) is used as an observation-based estimate of the stratification index. The same horizontal area-weighting treatment as for the ESMs is applied to the three-dimensional ⍴ field of WOA13 leading to a finer verticallyresolved (102 levels) one-dimensional vertical profile. ⍴ anomaly profiles comparing the ESMs and WOA13 are computed by vertically interpolating the highresolution WOA13 ⍴ profile to each coarsely-resolved model levels. The standard deviation of the WOA13 climatological monthly mean ⍴ is used as a proxy for the uncertainty around the climatological mean, as such uncertainty is not provided in the WOA13 database 26 . Standard statistical formulas 56 for uncertainty propagation are applied for the three-to-one-dimensional reduction. The SI standard deviation (σ SI ) of the observational constraint is calculated from the SI formula: where σ ρ z 0 and σ ρ z i are the WOA13 standard deviations.

Code availability
The MATLAB environment was used for statistical processing, model analyses and figure creation. The Gibbs-SeaWater (GSW) Oceanographic Toolbox has been used to convert model sea potential temperature and practical salinity to conservative temperature and absolute salinity, and to calculate in situ density (http://www.teos-10.org/software.htm).