Decoupling of δ¹⁸O from surface temperature in Antarctica in an ensemble of historical simulations

Stable water isotopes recorded in Antarctic ice cores have traditionally been used to infer past surface air temperatures (SATs). During the historical period (1850 onward), observational data and good-quality ice core records overlap, yielding an opportunity to investigate key relationships between ice core stable water isotope (δ¹⁸O) measurements and the Antarctic climate. We present a new ensemble of climate model simulations covering 1851–2004 using the UK Met Office HadCM3 general circulation model equipped with stable water isotopes. Our ensemble captures observed historical SAT and precipitation trends and weak δ¹⁸O trends. The weak δ¹⁸O trends mean there is no significant relationship between SAT and δ¹⁸O over one-third of Antarctica, and also half of our considered ice core sites, though relationships are stronger when using regional averages. The strongest regional relationships occur in the West Antarctic Ice Sheet (WAIS) region. This decoupling between SAT and δ¹⁸O occurs primarily because of the impact of autumnal sea ice loss during the simulated warming. The warming and sea ice loss are associated with (i) changes in near-coastal air mass intrusions (synoptic effects) induced by changes in the large-scale circulation and/or sea ice; (ii) direct sea-ice-driven changes in moisture pathways (especially lengths) to Antarctica; and (iii) precipitation seasonality changes, again mostly driven by sea ice changes. Consequently, when reconstructing temperatures over these timescales, changes in sea ice need to be considered, both to determine the most appropriate SAT and δ¹⁸O relationship and to understand how uncertainties affect the inference of past temperature from ice core δ¹⁸O measurements.

1 Introduction

Strong visible signs of Antarctic response to climate change have recently emerged. While a new sea ice cover minimum was recorded in February 2022 (Raphael and Handcock, 2022; Turner et al., 2022), with an extent of 1.97 million km², this record was broken the following year with sea ice extent falling to 1.91 million km² on 13 February 2023, associated with strong westerly winds and a 1.5 °C positive anomaly for Antarctic Peninsula air temperatures. The collapse of Antarctic ice shelves has similarly increased in frequency (Graham et al., 2022; Milillo et al., 2022; Wille et al., 2022). The consequent weakening of the buttressing force from the sea-ice-free areas and ice shelf collapse acts to accelerate Antarctic ice loss. Thus, the warming of Antarctica will have significant consequences for the global and regional mean sea level (Edwards et al., 2021; Seroussi et al., 2020; Parsons et al., 2020; Garbe et al., 2020) alongside consequences for Antarctic life and its environs (Golledge et al., 2019; Post et al., 2019), which require thought about adaptation (IPCC, 2022).

The relatively short satellite record (only since 1979) and the sparsity of in situ observational data from Antarctica mean that reconstruction of past temperature change is important for understanding natural variability and hence for our ability to detect anthropogenic climate change in Antarctica (Turner et al., 2004; Casado et al., 2023). Our understanding of pre-industrial climate change and its variability is mostly based on the reconstruction of temperature from proxy data (PAGES 2k Consortium, 2019). In Antarctica, stable water isotopes are the measurement most commonly used to reconstruct past surface air temperatures (SATs). This type of reconstruction is generally based on an empirical relationship between present-day (PD) surface snow water isotopes and surface air temperature (Dansgaard, 1953). The relationship between SAT and the ratio of heavy-water to light-water isotopes (expressed as δ¹⁸O) from Antarctic surface snow is usually assumed to be linear. When this linear relationship is used to estimate past SAT values from ice core water isotope measurements, it is sometimes referred as the “isotopic palaeothermometer” (Lorius et al., 1969; Masson et al., 2000).

The isotopic palaeothermometer approach has been successfully applied to deep ice cores to reconstruct past temperatures on long timescales (Jouzel et al., 2007; Lambert et al., 2008; Masson-Delmotte et al., 2010; Wolff et al., 2010). Casado et al. (2023) recently reconstructed the past 1000-year temperature record using an isotope-to-temperature conversion. This ice-core-based record was then used to show that the simulated temperature variability from the atmospheric general circulation models (AGCMs) run in the frame of the Coupled Model Intercomparison Project (CMIP) phases 5 (Taylor et al., 2012) and 6 (Eyring et al., 2016) is too low. The isotopic palaeothermometer relationship has been shown to vary spatially over the Antarctic continent (e.g. Sime et al., 2008, 2009a). Goosse et al. (2012), PAGES 2k-PMIP3 group (2015), and Neukom et al. (2018) show that noise and the spatial coverage of δ¹⁸O and other proxy data affect our understanding of these variations, while Goursaud et al. (2018, 2019) and others show smaller, or less reliable, changes in δ¹⁸O for a given temperature change in coastal regions (Isaksson and Karlén, 1994; Sime et al., 2008, 2009a; Abram et al., 2013; Thomas et al., 2013; Goursaud et al., 2017). Data to investigate SAT–δ¹⁸O relationships are sparse (Masson-Delmotte et al., 2008; Landais et al., 2017); nevertheless, the PAGES 2k Network Antarctica2k (A2k) helped to address this question of geographical variations in the SAT–δ¹⁸O relationship by defining regions and compiling the available δ¹⁸O data from ice cores (Stenni et al., 2017).

The geographical variability in the isotopic palaeothermometer is due to controls on δ¹⁸O rather than related to SAT. These other controls include changes related to atmospheric dynamics, such as changes in the synoptic and seasonal nature of precipitation (van Ommen and Morgan, 1997; Krinner and Werner, 2003; Jouzel et al., 2003; Sime et al., 2008; Servettaz et al., 2023b) and air mass sources (Landais et al., 2021); various impacts from changes in Antarctic ice sheet morphology (Holloway et al., 2016; Werner et al., 2018; Buizert et al., 2021; Goursaud et al., 2021); and sea ice variability (Holloway et al., 2018; Cauquoin et al., 2023). The stability of the SAT–δ¹⁸O relationship has thus been of much interest for more than 2 decades (Jouzel et al., 2003). Following Sime et al. (2008) and Sime et al. (2009a), the importance of changes in synoptic events in the context of the anthropogenic warming was recently explored by Wille et al. (2019) and Pohl et al. (2021), who confirm that the impact of synoptic changes on the SAT–δ¹⁸O relationship can be important for the past palaeothermometer during warm climates (Dalaiden et al., 2020).

AGCMs equipped with stable water isotopes are a key tool to investigate the climate processes driving temporal variability in the palaeothermometer relationship (e.g. Werner et al., 2001; Sime et al., 2008; Werner et al., 2018). For instance, AGCM isotopic studies have focused on the effects of external forcing on the SAT–δ¹⁸O relationship, including elevation and greenhouse gases across a range of timescales (e.g. Sime et al., 2009b; Werner et al., 2018; Goursaud et al., 2021). A major result is that, for differing timescales and driving mechanisms, different SAT–δ¹⁸O relationships can be obtained. This emphasises the importance of investigating the impact of atmospheric dynamical drivers, particularly changes in sea ice and precipitation seasonality during past warm periods in Antarctica (Sime et al., 2008, 2009b; Holloway et al., 2016). Only one AGM study has investigated the signature of stable water isotopes in historical simulations (Yoshimura et al., 2008). While the applicability of the Antarctic palaeothermometer relationship has been investigated across various timescales, it has not yet been thoroughly investigated over the historical period. Furthermore, it has not been investigated using transient historical (1851–2004) simulations.

Here, we run an ensemble of transient historical (1851–2004) simulations, using the stable-water-isotope-enabled coupled general circulation model, HadCM3. This ensemble provides a benchmark of historical precipitated stable water isotopes covering the whole continent and allows the SAT–δ¹⁸O relationship over the historical period to be investigated. Firstly, we examine trends in SAT, precipitation, and δ¹⁸O and compare these against observed trends. We then examine SAT–δ¹⁸O relationships, including regional patterns, and the question of model dependence. Finally, in order to understand the drivers of δ¹⁸O change, we perform a detrended composite analysis for cold- and warm-year ensembles and quantify the impact of seasonality changes in precipitation and δ¹⁸O. This analysis provides an understanding of SAT–δ¹⁸O variability and, in particular, the role of sea ice change in this relationship during the historical period.

2 Materials and methods

2.1 Model and simulations

Here, we use the Hadley Centre atmosphere–ocean general circulation model (AOGCM), HadCM3, to run six transient historical simulations, HadCM3 is a version of the coupled atmosphere–ocean UK Met Office climate model (Pope et al., 2000; Gordon et al., 2000), which means that sea ice is prognostic.. The model is equipped with stable water isotopes (Tindall et al., 2009). Its horizontal resolution is 3.75°×2.5°, and there are 19 vertical levels (Pope et al., 2000; Gordon et al., 2000; Tindall et al., 2009).

The setup of the historical simulations is described in Schurer et al. (2014) and follows the recommendations of the third Paleoclimate Modelling Intercomparison Project (PMIP3; Schmidt et al., 2011). Each simulation is forced with time-varying orbital, solar, volcanic, land-use, and well-mixed greenhouse gas forcing. As above, sea ice is not prescribed but rather calculated by the model. Changes in orbital parameters were calculated following Berger (1978). Volcanic forcing is that described in Crowley et al. (2008). The solar forcing follows Shapiro et al. (2011). Changes in CO₂, N₂O, and CH₄ were set following the PMIP3 standard (Schmidt et al., 2011). Changes in the abundances of six halocarbons were prescribed following Tett et al. (2007). Changes in land cover were prescribed by reclassifying the global land cover reconstruction developed by Pongratz et al. (2008). Each of our simulations were only altered by starting each simulation 1 year apart.

We analyse HadCM3 surface air temperature (SAT; °C), precipitation (P; mm per month), precipitation-weighted δ¹⁸O (δ¹⁸O), and sea ice extent (defined as the region of ice-covered ocean where the sea ice concentration is >15 %). HadCM3 provides a reasonable representation of Antarctic climate and δ¹⁸O (Appendix A and Turner et al., 2006; Tindall et al., 2009; Holloway et al., 2016).

2.2 Data and methods

We perform a model–data comparison using the Stenni et al. (2017) ice core data compiled by the PAGES A2k project by binning our model output, including δ¹⁸O, into 5-year-equivalent averages and compute anomalies relative to the 1960–1990 mean. In order to investigate the historical mean climate state and variability, we compute ensemble mean values covering the period 1851–2004 using monthly outputs. Trends over the whole period and for the last 50 years are calculated using linear regressions. Where we regress climate variables against δ¹⁸O, the linear regressions are computed using the stacked individual ensemble members rather than using the ensemble mean. This approach ensures that the ensemble variability is included in our linear regression statistics and increases the number of points on which the regressions are processed. Gradients from the linear regressions are provided with a plus/minus standard error. Results from linear relationships are stated only where they are significant, using a p-value ≤ 0.05.

Our historical SAT–δ¹⁸O linear relationship at the regional scale is compared with the regional slopes and correlation coefficients that we computed from the AGCM ECHAM6-wiso equipped with water stable isotopes (Cauquoin et al., 2019). The water stable module of this last generation of the model ECHAM was updated compared to its predecessor, especially (i) the supersaturation parameters; (ii) the kinetic fractionation at the evaporation over oceans, now assumed to be independent of the wind speed in order to better represent the d-excess versus deuterium relationship from the Antarctic Snow reported by Masson-Delmotte et al. (2008); and finally (iii) the sublimation processes now accounting for the isotopic content of snow over sea ice. Here, we use a simulation run at a T127L95 resolution ($\sim \mathrm{0.9}\mathit{°}\times \mathrm{0.9}\mathit{°}$ horizontal resolution and 95 vertical levels) and nudged towards the ERA5 reanalyses (Hersbach et al., 2020) over the period 1979–2022 (Cauquoin and Werner, 2021).

Composites are used to interpret our results. Warm and cold (versus mean) annual composites results are defined using detrended annual area-weighted SAT. The years of the ensemble mean with SAT below (above) the mean minus (plus) 1 standard deviation constitute a cold (warm) ensemble.

To examine the impact of changing seasonality over the historical period, we isolate the impact of precipitation and δ¹⁸O seasonal changes, recorded in the ensemble mean, on the precipitation-weighted δ¹⁸O between the first 50 years of the simulation and the last 50 years of the simulations (see Liu and Battisti, 2015; Holloway et al., 2016; Sime et al., 2019). This is calculated as

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{{\mathrm{\Delta }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{seas}}={\displaystyle \frac{\sum _j{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_j^{\mathrm{recent}}\times P_j}{\sum _j{P}_j}}-{\displaystyle \frac{\sum _j{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_j\times P_j}{\sum _j{P}_j}},\text{(2)}& {\displaystyle }& {\displaystyle }{P}_{\mathrm{seas}}={\displaystyle \frac{\sum _j{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_j\times P_j^{\mathrm{recent}}}{\sum _j{P}_j^{\mathrm{recent}}}}-{\displaystyle \frac{\sum _j{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_j\times P_j}{\sum _j{P}_j}}.\end{array}$$

The summations, with index j, are over the 12 months of the year. Variables with superscript “recent” indicate that they were extracted for the last 50 years of the simulation, whereas the variables without superscript indicate that the variables were extracted for the first 50 years of the simulation.

3 Results

This section uses these model data and methods to examine trends in Antarctic SAT, precipitation, sea ice, and δ¹⁸O, including at the continental and regional scale; relationships between temperature versus δ¹⁸O, including their stability and model dependency; and, finally, the drivers of δ¹⁸O changes.

3.1 Trends in Antarctic SAT, precipitation, sea ice, and δ¹⁸O

We analyse our simulations against available observations and reanalysis data, e.g. the Climate Forecast System Reanalysis (CFSR) and the National Centers for Environmental Prediction reanalyses 2 (NCEP-2), at the Antarctic scale before evaluating regional-scale climate changes. This is followed by a comparison of simulated δ¹⁸O changes against the Stenni et al. (2017) ice core dataset. We note that all trends outlined below are similar regardless of whether we use the full historical period or the last 50 years of the ensemble simulations. Nevertheless, to permit the most direct comparison where possible, we have matched our calculations to the periods used by other authors.

3.1.1 Continental trends in climate

Our simulated SAT trend over Antarctica is 0.12±0.02 °C per decade over the last 50 years (Fig. 1). This is consistent with observations of 0.12±0.07 °C per decade over 1957–2006 (Steig et al., 2009) and 0.11±0.08 °C per decade over 1959–2012 (Nicolas and Bromwich, 2014). Our trend is lower than the 0.22±0.04 °C trend of Casado et al. (2023); however, this is itself dependent on their isotope-to-temperature reconstruction method. Over East Antarctica, our simulated historical SAT trend is weaker (0.10±0.02 °C per decade; r=0.59) than over West Antarctica (0.15±0.02 °C per decade; r=0.75). This compares with 0.10±0.07 °C per decade for the east and 0.17±0.06 °C per decade for the west calculated by Steig et al. (2009).

Figure 1 Antarctic trends over the periods 1851–2004 and 1955–2004. Each column is associated to a climate variable: (a, e) surface air temperature (SAT in °C, in red), (b, f) precipitation (P in mm per month, in blue), (c, g) δ¹⁸O (in ‰, in green), and (d, h) sea ice concentration (SIC in %, in grey). Panels (a)–(d) display the time series of Antarctic surface-weighted averages over the period 1851–2004. Solid coloured lines represent the annual average of the ensemble mean, the coloured surfaces represent the annual standard deviations of the ensemble mean, the solid grey lines represent the simulations, and the dashed lines represent the linear regressions. The slopes (a) plus or minus the standard errors and the correlation coefficients (r) are given in the top right of the figures, with the first row corresponding to the whole historical period 1851–2004 and the second row, in bold, to the last 50 years 1955–2004. Panels (e)–(h) display the maps of anomalies over the period 1955–2004 against 1854–1904.

Our simulated historical precipitation trend over the last 50 years is 3.1 mm yr⁻¹ per decade (Fig. 1), which lies between the admittedly wide Bromwich et al. (2011) reanalysis-based equivalent values of 0.4±1.8 mm yr⁻¹ per decade from the Climate Forecast System Reanalysis (CFSR), and 7.1±1.5 mm yr⁻¹ per decade from the National Centers for Environmental Prediction reanalysis 2 (NCEP-2) over 1979–2009. The HadCM3 results are also in agreement with Dalaiden et al. (2020), who show an increase in both West and East Antarctic Ice Sheet (WAIS and EAIS)-simulated precipitation.

Our HadCM3-simulated historical Antarctic September sea ice decrease is $-\mathrm{0.20}\pm \mathrm{0.01}\times {\mathrm{10}}^{\mathrm{6}}$ km² per decade over the period 1851–2004 ($r=-\mathrm{0.69}$) and $-\mathrm{0.40}\pm \mathrm{0.06}\times {\mathrm{10}}^{\mathrm{6}}$ km² per decade over the period 1954–2004 (r=0.67). This is consistent with the recent results of Shu et al. (2020), who calculated Antarctic September sea ice trends of −0.45 and $-\mathrm{0.43}\times {\mathrm{10}}^{\mathrm{6}}$ km² per decade during the period 1979–2005 from the Coupled Model Intercomparison Project (CMIP) 5 and 6 model results respectively. Whilst neither the Shu et al. (2020) nor our HadCM3 sea ice trends match the observed slope of $+\mathrm{0.10}\times {\mathrm{10}}^{\mathrm{6}}$ km² per decade (p=0.16) (Shu et al., 2020), they do agree with each other. It is also worth noting that, since 2016, Antarctic sea ice has begun losing significant area during both summer and winter.

The simulated historical changes in Antarctic-wide SAT and precipitation and sea ice thus seem in reasonable agreement with other model results and possibly also observations (though this is less clear). Trends, however, vary between regions.

3.1.2 Regional-scale trends

We now look at the simulated climate, using the Antarctic regions defined by the PAGES A2k community (see Fig. 2, regions defined by Stenni et al., 2017). Every Antarctic region shows a simulated increase in SAT over the historical period (1851–2004; Fig. 2), but these vary across the continent. Warming trends are strongest for the Peninsula and Dronning Maud Land regions (0.11 °C per decade and r≥0.8), closely followed by the WAIS and the Weddell coast regions (0.08 °C per decade and r≥0.8). The Plateau and the Indian coast show weaker warming trends of 0.04 and 0.05 °C per decade respectively (r=0.6 for both). These historical trends approximately match the Turner et al. (2020) results derived from station data.

Figure 2 Regional trends over the periods 1851–2005 and 1955–2005. Time series of surface-weighted averages of (a) surface air temperature (SAT in °C, in red), (b) precipitation (P in mm per month, in blue), (c) δ¹⁸O (in ‰, in green), and (d) sea ice concentration (SIC in %, in grey) over the period 1851–2005 for the different Antarctic regions as defined by Stenni et al. (2017). Solid coloured lines represent the annual average of the ensemble mean, the coloured surfaces represent the annual standard deviations of the ensemble mean, the solid grey lines represent the simulations, and the dashed lines represent the linear regressions. The slopes (a) and the correlation coefficients (r) are given in the top right of the figures. The first row reports the regional slope and correlation coefficient over the whole historical period, while the second row reports the regional slope and correlation coefficient over the last simulated years. Shaded grey rows correspond to non-significant relationships (p-value > 0.05).

Our simulated SAT trends resemble the Stenni et al. (2017) regional warming trends. Stenni et al. (2017) found cooling trends for the Plateau, the Weddell coast, and Victoria Land (from −0.13 to −0.05 °C per decade), where we found a small warming. They found a larger trend for the Peninsula (from 0.2 to 0.29 °C per decade), though, again, we point out that the Stenni et al. (2017) values, like the Casado et al. (2023) results, are partly dependent on isotope-to-temperature reconstruction methods. At the scale of station locations, Jones et al. (2019) also show the highest trends for the Peninsula.

Similarly to SAT, HadCM3 shows an increase in precipitation for all the Antarctic regions, both for the full historical period and for the last 50 years. Similarly to temperature, the Peninsula features the strongest trend of 7.8 mm yr⁻¹ per decade over the historical period, whereas the Plateau and the Ross sections display weaker trends of 0.5 and 0.82 mm yr⁻¹ per decade over the same period. Interestingly, Thomas et al. (2017) and Medley and Thomas (2019) found similar results. Whilst these results are from ice cores, they are not dependent on the interpretation of isotopes; instead, they mostly use profiles of density and layer counting from relevant age markers (e.g. chemical species, radio isotopes, biologically compounds).

For sea ice, HadCM3 simulates a sea ice decrease around all sectors of Antarctica, except for the Weddell sector (Table B1). Trends are largest in the Indian sector ($-\mathrm{49}\pm \mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ km² per decade; $r=-\mathrm{0.67}$) and smallest in the Pacific sector ($-\mathrm{25}\pm \mathrm{4}\times {\mathrm{10}}^{\mathrm{3}}$ km² per decade; $r=-\mathrm{0.44}$). Whilst these results are not compatible with pre-2016 satellite observations, they are consistent with other climate simulations (Shu et al., 2020).

Despite the simulated increases in SAT and precipitation, δ¹⁸O simulated by HadCM3 shows a very weak trend of 0.04±0.003 ‰ per decade (r=0.21) over the last 50 years. Interestingly, Casado et al. (2023) provide a higher trend from 1950–2005 of 0.11±0.02 ‰ per decade, based on ice core data. Different reasons, which we are not able to elucidate so far, could explain that mismatch, inter alia, (i) a model discrepancy to resolve processes, (ii) the model resolution, (iii) the geographical distribution of the ice core locations, and (iv) the different methods for the SAT–δ¹⁸O calibration. Section 5 of this paper focuses on investigating and explaining the δ¹⁸O trends. Before this, we provide a brief overview of the regional picture.

At the regional scale, over the historical period, trends are small (Fig. 2). In terms of linear relationship, it is null for Victoria Land, while the gradient is the highest for the Weddell coast with a trend of 0.05 ‰ per decade (r=0.39) and the correlation coefficient is the highest (e.g. the strongest linear relationship) for the Peninsula with a trend of 0.04 ‰ per decade (r=0.57). Over the last 50 years, only three regions, the Indian, Weddell and Dronning Maud Land coastal regions, keep on displaying significant δ¹⁸O trends, which double or more compared to the historical period, with gradients of 0.08 ‰ per decade, 0.08 ‰ per decade, and 0.14 ‰ per decade respectively. Stenni et al. (2017) compiled δ¹⁸O trend statistics based on ice core anomalies using unweighted composites over the period 1900–2000, based on 5-year bins. They found only three regions with significant trends, which are the Indian coast, the Peninsula, and Dronning Maud Land, with gradients in the range of our results for the Indian coast and higher gradients for the Peninsula and Dronning Maud Land (mean trends of 0.15 ‰ per decade and 0.11 ‰ per decade respectively). Comparatively, Casado et al. (2023) calculated trends over windows varying between 35 and 60 years and using a persistence method. They found gradients with the same range of values, from 0.09 ‰ for the Indian coast to 0.19 ‰ for the Weddell coast, while they found significant relationships where we do not for time windows varying from 40 to 65 years. Note that, for most of the regions, the significance of our simulated relationships disappears for time windows shorter than 75 years (Appendix D). This could be explained either by the simulated anthropogenic variability being too low, as suggested by Casado et al. (2023), or by a change in the drivers on δ¹⁸O. The disparities between our results and the previous studies could be explained by the different time windows, the different methodologies, the lack of ice core data to make representative regional reconstructions, or a model discrepancy. While Casado et al. (2023) carefully investigated the impact of the data stack method and the time window on the δ¹⁸O reported trends, we suggest that an extended study could compare the statistical and dynamical methods on both ice core data and water-stable-isotope-enabled AGCM outputs to complete the analysis.

3.2 Temperature versus δ¹⁸O relationships

Given that much of ice core science is underpinned by the relationship between temperature and δ¹⁸O (Jouzel et al., 2003), and having discussed simulated SAT, precipitation, sea ice, and δ¹⁸O trends, we now investigate temperature versus δ¹⁸O relationships. Sime et al. (2008) demonstrate that there is no clear relationship between the spatial versus temporal SAT–δ¹⁸O gradients across Antarctica. We therefore focus on comparing the last 50 years of our simulations against the whole historical period. To enable a consideration of model dependency, we also compare our historical ensemble against a nudged ECHAM6-wiso simulation (Table 1).

Table 1 Historical SAT–δ¹⁸O relationships at the regional scale. Slope (in ‰ °C⁻¹) plus or minus the standard error and the correlation coefficient (in brackets) of the surface-weighed average of surface air temperature against the surface-weighed average of δ¹⁸O for the Antarctic regions as defined in the PAGES Antarctica2k project (Stenni et al., 2017): the Plateau, the Indian coast, the Weddell coast, the Peninsula, the WAIS, Victoria Land, and Dronning Maud Land. These are simulated by the ECHAM6-wiso model (over the period 1979–2022, 44 points, “ECHAM6-wiso”), by HadCM3 over the last 50 years (1955–2004, 50 points, “last 50 years of HadCM3”), and by HadCM3 over the whole historical simulated period (1851–2004, 154 points, “historical HadCM3”) using the ensemble mean of the six simulations (see Materials and methods). All the relationships are significant (p-values < 0.05).

3.2.1 Antarctic-wide and regional-scale

The simulated SAT–δ¹⁸O relationship, calculated using annual means on each grid point, is statistically significant over 66 % of Antarctica (Fig. 3). For the continent as a whole, we simulate a mean Antarctic SAT–δ¹⁸O gradient of 0.57±0.06 ‰ °C⁻¹ (r=0.62) over the historical period, increasing to 0.67±0.13 (r=0.60) over the last 50 years of the simulation. The comparable numbers for an Antarctic-wide SAT–δ¹⁸O relationship in Casado et al. (2023) are 0.49 ‰ °C⁻¹ to 0.69 ‰ °C⁻¹. In our simulations, non-significant SAT–δ¹⁸O relationships occur in eastern Antarctica between 40 and 100° E; in all regions between 140 and 220° E covering the Wilkes coast, Victoria Land, and some parts of Queen Mary Land; and on the coast of Dronning Maud Land with some areas at 350–360° E joining the South Pole. Non-significant relationships were also reported in observations and model outputs. For instance, Goursaud et al. (2018) report no SAT–δ¹⁸O relationship at the annual scale over the coast of Dronning Maud Land, Victoria Land, some of the Indian coast, and the Peninsula. Results showing the absence of an SAT–δ¹⁸O relationship derived from firn/ice cores were also published (e.g. Goursaud et al., 2019; Bertler et al., 2011; Vega et al., 2016; Goursaud et al., 2017). More recently, Casado et al. (2023) also found no significant relationships, though possibly because of a lack of data, for some regions. For them, these regions were the Indian and Weddell coasts and Victoria Land (see section 3.3 for a discussion of why these regions do not show statistically significant temperature versus δ¹⁸O relationships.) In contrast, the Peninsula, part of the WAIS coast, and some parts of the Plateau show much stronger SAT–δ¹⁸O relationships. Indeed, some coastal areas are associated with the highest correlation coefficients, ranging between 0.15 and 0.45. The highest SAT–δ¹⁸O gradients in the Plateau region and south of the Filchner Ice Shelf can exceed 0.75 ‰ °C⁻¹. Casado et al. (2023) also found the largest SAT–δ¹⁸O gradients in similar regions.

Figure 3 Historical SAT–δ¹⁸O pattern. The SAT–δ¹⁸O slopes (a) and correlation coefficients (b) simulated by HadCM3 over the historical period at the interannual scale at each grid point. Gradients (in  ‰ °C⁻¹, A) and correlation coefficients (r) of the SAT–δ¹⁸O relationships at the interannual scale (151 points) at each Antarctic grid point and based on the stacked simulations. Regions with hashed black lines indicate the presence of no significant relationships (p-value ≥ 0.05).

This Antarctic-wide picture of geographical variability in the temperature versus δ¹⁸O relationship is reasonably consistent with previous studies covering parts of the historical period (Sime et al., 2009a; Stenni et al., 2017; Goursaud et al., 2019) and measurements from coastal firn cores (Isaksson and Karlén, 1994; Abram et al., 2013; Thomas et al., 2013; Goursaud et al., 2017). Interestingly, Guan et al. (2016, 2020) associate similar results (low or negative SAT–δ¹⁸O relationships) with higher source temperatures. Here, we suggest this is mainly driven by sea ice retreat (see Sect. 3.3).

3.2.2 Stability over the historical period and model dependency

Results from HadCM3 are similar for both the last 50 years and the whole historical period (Table 1): SAT–δ¹⁸O gradients vary between 0.3 ‰ °C⁻¹ and 0.7 ‰ °C⁻¹. The average difference is 0.07 ‰ °C⁻¹. The only region with a statistically different result is Dronning Maud Land, with SAT–δ¹⁸O gradients of 0.76±0.12 ‰ °C⁻¹ and 0.49±0.05 ‰ °C⁻¹ over the last 50 years and the whole historical period respectively. Thus, for most of the continent, our HadCM3 results over the last 50 years and the whole historical period appear equivalent.

Figure 4 Composite maps for cold (a–d) and warm (e–h) years. Maps of Antarctic (a, e) δ¹⁸O (in ‰), (b, f) surface air temperature (SAT in °C), (c, g) precipitation (P in mm per month), and (d, h) sea ice concentration (SIC in %) for years with surface air temperatures below (a–d) and above (e–h) 2 standard deviations in the ensemble mean over the period 1851–2004. Maps use the detrended ensemble mean (see Materials and methods). There are 8 years (out of the 155 simulated years) with SAT anomalies out of the plus or minus 2 standard deviations: 4 are above that range, and 4 are below that range. Panels (c) and (g) also show the wind field. Blue (red) arrows indicate southward (northward) winds, and regions of the southward winds are delimited using red contours.

Interestingly, the ECHAM6-wiso SAT–δ¹⁸O gradients calculated here are on average 2 times lower than those computed using the ECHAM5-wiso simulations nudged to ERA-interim over the period 1979–2013 and published in Stenni et al. (2017). Thus, the improvements to the ECHAM6 coding seem to bring HadCM3 and ECHAM into alignment (Table 1). We consider these new ECHAM6-wiso and HadCM3 values to likely be more accurate. We can note the small differences between HadCM3 and ECHAM6-wiso: ECHAM6-wiso simulates slightly stronger relationships with a mean correlation coefficient difference of 0.04, while gradients tend to be slightly higher in HadCM3 with a gradient difference of 0.13 ‰ °C⁻¹. The only notable differences are for Dronning Maud Land and the Indian coast, with stronger relationships and higher gradients simulated by HadCM3. Thus, whilst it is unclear whether the nudging of ECHAM6-wiso towards ERA5 reanalysis, the model resolution, the model physics, or the difference in sea ice behaviours is the main reason for these discrepancies, it is clear that simulated temperature versus δ¹⁸O relationships have low but significant uncertainties. These need to be considered, both regionally and for the most relevant climate state, before undertaking any inferences of past temperatures using isotopes measured in ice cores.

3.3 Drivers of δ¹⁸O changes

We use two approaches to investigate the mechanisms driving simulated δ¹⁸O changes. Firstly, we separate and compare extremely warm and cold years for both annual (Fig. 4, Table C1) and seasonal (Fig. 5) data by generating annual and seasonal composites with mean annual Antarctic SAT anomalies greater than plus or minus 2 standard deviations from the mean respectively. Secondly, we isolate the impact of changing precipitation seasonality on δ¹⁸O, showing simple month values (Fig. 6) and also following the decomposition method used in Liu and Battisti (2015), Holloway et al. (2016), and Sime et al. (2019) (Fig. 7).

Figure 5 Seasonal warm–cold changes. Differences between the warm and cold ensemble means for winter (i.e. June to August, a–c) and summer (i.e. from December to February, d–f) surface air temperature (a, d, SAT in °C), precipitation (b, f, P in mm per month), and precipitated δ¹⁸O (c, f, δ¹⁸O in ‰).

Figure 6 Impact of climate seasonal changes. Seasonal differences in the precipitation (in %), the precipitated δ¹⁸O (in ‰), and the sea ice concentration (in %) between the first 50 simulated years and the last 50 simulated years. Light lines correspond to the member simulations, and the dark line corresponds to the ensemble mean.

Figure 7 The impact of seasonal changes on precipitation and on mean annual δ¹⁸O. δ¹⁸O changes (Δδ¹⁸O in ‰) due to changes in the seasonal cycle of precipitation (a) and due to changes in the seasonal cycle of δ¹⁸O (b) between the last 50 simulated years and the first 50 simulated years of the HadCM3 historical mean ensemble. See Materials and methods for details of the decomposition.

As expected, the spatial patterns of SAT and sea ice anomalies tend to vary together, with the pattern being approximately mirrored between cold and warm composites (Fig. 4, top and bottom panels respectively). Whilst fully isolating the drivers of δ¹⁸O is tricky, Figs. 4 to 7 suggest that the primary mechanism driving continental-scale SAT–δ¹⁸O decoupling in HadCM3 is the simulated loss of sea ice over the historical period (Fig. 5d and h).

The September average sea ice area across the warm composite is 5.8×10⁶ km² less than in the cold composite. Given that this reduction occurs primarily during winter (Fig. 5c; there is almost no summertime sea ice around Antarctica), warmer years tend to receive relatively more precipitation during winter months compared to cold years, partially offsetting the warming signal in δ¹⁸O. This can be seen in Fig. 5, displaying seasonal anomalies (for winter, e.g. from June to August, and for summer, e.g. from December to February) in precipitation, δ¹⁸O, and sea ice between the warm and cold composites: the largest (smallest) precipitation and δ¹⁸O anomalies occur during the winter (summer) months. Precipitation anomalies peak in autumn and winter, whilst δ¹⁸O anomalies peak in winter and spring (Fig. 6), the latter coincident with the annual maximum sea ice extent and largest sea ice area anomalies. The relative increase in winter precipitation during warm years acts to reduce δ¹⁸O across Antarctica, compared to if the seasonality of precipitation remained unchanged. This is perhaps seen most clearly in Fig. 7, where Fig. 7a is predominantly blue, which indicates that precipitation seasonality changes are acting to decrease δ¹⁸O. The effect of changing seasonality is particularly large in the Indian, Dronning Maud Land, and Victoria Land (through the Wilkes Land) sectors, which are prone to air mass intrusions (Figs. 5c and 7a).

Although the reduction in sea ice area simulated in the warm composite and throughout the historical ensemble increases the proportion of winter precipitation, negatively influencing δ¹⁸O, less sea ice also shortens the distance between evaporation source and precipitation site, which has an opposing positive influence on δ¹⁸O. This effect is evident as negative sea ice anomalies adjacent to coastal regions with large positive δ¹⁸O anomalies in Fig. 7b and in Fig. 5c and f. Sea ice loss may also allow locally sourced precipitation to penetrate further inland into the Antarctic interior, which usually receives precipitation sourced from lower latitudes (Gao et al., 2024), and promote an increase in high-intensity precipitation events during cold seasons (e.g. Schlosser et al., 2004). Consequently, sea ice loss during the historical period leads to competing influences on δ¹⁸O and considerable spatial variability in SAT–δ¹⁸O relationships (Fig. 3a), with no significant relationship in several regions. Changes in the precipitation seasonality (Fig. 7a) reduce mean Antarctic δ¹⁸O by 0.16±0.19 ‰ and for the lowest changes by 0.11  ‰ over the Plateau and no changes over the Peninsula. Changes in the seasonal cycle of δ¹⁸O (Fig. 7b) are more spatially variable and largely sea-ice-driven, with a depletion up to −0.45 ‰ over the Plateau and Victoria Land and an enrichment from 0.1 ‰ to 1.0 ‰ in coastal regions.

These historical simulations indicate that Antarctic δ¹⁸O is highly sensitive to patterns of sea ice change, which influence atmospheric dynamics, air mass pathways and lengths, and water isotope evaporation and condensation temperatures. Our results are particularly sensitive to autumnal sea ice changes: the largest simulated reduction in sea ice occurs in autumn, coincident with the largest changes in δ¹⁸O (Fig. 5c). The dynamic processes behind the δ¹⁸O changes induced by sea ice extent are complex and multiple. Although the Southern Annular Mode, the leading mode of the atmospheric variability in the Southern Hemisphere, might explain some of these δ¹⁸O-simulated changes (Appendix G), a more comprehensive study might investigate the impact of the atmospheric circulation changes.

4 Conclusions

Results from six transient simulations over the historical period allow us to examine Antarctic precipitation-weighted δ¹⁸O and its relationship with Antarctic SAT. The ensemble features a rise in mean Antarctic SAT of 0.12±0.02 °C per decade over the last 50 years, consistent with previous studies (Steig et al., 2009; Nicolas and Bromwich, 2014). In agreement with observations, the associated simulated trend in the water isotope δ¹⁸O is weak ($+\mathrm{0.04}\pm \mathrm{0.03}$ ‰ per decade). Unlike SAT, the δ¹⁸O trend is weaker during the last 50 years compared to the complete historical period. This implies that, over the last 50 years of our ensemble, δ¹⁸O is influenced by processes other than SAT-related condensation temperature; i.e. ice and atmospheric processes oppose purely SAT and condensation temperature controls on δ¹⁸O.

Although the consequences of these competing effects vary spatially, they result in non-statistically significant relationships between SAT and δ¹⁸O for approximately one-third of the continent over the historical ensemble. Non-significant SAT–δ¹⁸O relationships occur across much of Antarctica between 40 and 100° E: the Wilkes coast, Victoria Land, part of Queen Mary Land, the coast of Dronning Maud Land, and areas of the South Pole region. Interestingly, we find similar but slightly weaker SAT–δ¹⁸O correlations and slightly higher gradients compared to ERA5-nudged ECHAM6-wiso simulations at the regional scale.

We suggest three processes that lead to weak SAT–δ¹⁸O relationships during the historical ensemble. Firstly, historical changes in near-coastal air mass intrusions, induced by changes in the large-scale circulation and/or changes in synoptic events, have an impact on both the moisture source and the precipitation regime. Secondly, changes in sea ice concentration impact moisture pathways to Antarctica. Thirdly, changes in precipitation seasonality, largely driven by sea ice changes, bias the relative impact of cold- vs. warm-season precipitation on precipitation-weighted δ¹⁸O. Fyke et al. (2017) show similar spatial patterns in Antarctic precipitation during the pre-industrial (Fyke et al., 2017, Fig. 2c), driven by the large-scale circulation and acting to regulate atmospheric moisture transport and regional sea ice variability (Fyke et al., 2017; Marshall et al., 2017; Raphael et al., 2019). Changes in Antarctic sea ice (Holloway et al., 2016) and moisture source changes (Landais et al., 2021) have also been proposed as responsible for driving the Last Interglacial δ¹⁸O peak. Additionally, Sime et al. (2008) show that (sea-ice-related) changes in the seasonal and synoptic distribution of Antarctic precipitation were jointly responsible for a partial decoupling between SAT and δ¹⁸O. These results support the role of environmental changes, particularly sea ice, generating significant variability in the Antarctic SAT–δ¹⁸O relationship.

In conclusion, the results from this isotope-enabled historical ensemble permit investigation of the SAT–δ¹⁸O relationship during the historical period and the mechanisms driving its spatial variability. Whilst we recognise the limitations, in terms of spatial resolution and a simplistic sea ice sub-model, of our chosen model, HadCM3, we value the ability to perform a several-member ensemble of >100-year simulations using a fully coupled isotope-enabled model.

For future work, higher-resolution and more physically advanced models could be used, alongside new tracer methods (Gao et al., 2024). Firn models could also be used to examine the interesting question of post-deposition effects, such as the redistribution of snow by the wind (Libois et al., 2014), snow–vapour exchanges (Casado et al., 2016; Ritter et al., 2016), and snow metamorphism (Casado et al., 2021), on historical δ¹⁸O. Finally, more stable water isotope records from Antarctic ice and firn core data are more than needed to evaluate models and to lead model–data investigations of past climates, comparing SAT–δ¹⁸O relationships from different water-stable-isotope-enabled models, in line with the work of the Stable Water Isotope Intercomparison Group 2 (SWING) (IPCC, 2022).

Appendix A: HadCM3 evaluation of Antarctic surface climate and δ¹⁸O

A1 Method

Here, we check that HadCM3 provides a reasonable representation of the Antarctic surface climate and δ¹⁸O.

Surface air temperature (SAT) output data from HadCM3 are evaluated against the AntAWS dataset (Wang et al., 2023), a compilation of Antarctic observations from 267 automatic weather stations (AWSs) operational between some parts of the period from 1980 to 2021. Surface mass balance (SMB) model output, calculated within the model code as precipitation minus evaporation (wind-related processes are not accounted for by HadCM3), are similarly evaluated against AnSMB (Wang et al., 2021), the most recent quality-controlled published SMB compilation extracted from stakes, snow pits, ice cores, ultrasonic sounders, and ground-penetrating radar. Finally, simulated δ¹⁸O values are evaluated using the updated database compiled by Goursaud et al. (2018); this combines all available firn, ice core, surface snow, and precipitation observations of Antarctic δ¹⁸O.

We show maps and scatter plots (model versus observed values) for SAT, SMB, and δ¹⁸O. The comparison helps establish if the model underestimates the real spatial heterogeneity across Antarctica. Mean climatological values (20-year averages or more, averaged over the ensemble wherever possible) were calculated at each model grid point and directly compared to the most equivalent observational climatological value (see paragraph above). The comparison uses output from a closest grid point comparison method.

A2 Results

A2.1 SAT

The evaluation by Turner et al. (2006) of the HadCM3 Antarctic climate, especially including near-surface air temperatures, mean sea level pressures, and geopotential heights, shows a large warm bias in the Antarctic interior associated with a low-biased modelled orographic height (the highest model grid point elevations do not reach 4000 m a.s.l.). This finding remains fully consistent with the newer Wang et al. (2021) observation datasets (Fig. A1). The minimum climatological Antarctic Plateau SAT value is −37.2 °C (Fig. A1a), considerably warmer than the AntAWS minimum of −64.6 °C (Fig. A1c). In regions where the observational temperature is above −30 °C, the model values of SAT match the observations better, although there remains a slightly underestimate (warm bias) in West Antarctica (Fig. A1b and top right of Fig. A1c). Altogether, although the warm bias in the Antarctic interior contributes to weakening the linear regression between the HadCM3 simulations and the observations (correlation coefficient of 0.76), the Antarctic mean simulated SAT is surprisingly good: the Antarctic mean climatological SAT is $-\mathrm{25.1}\pm \mathrm{14.1}$ and $-\mathrm{25.0}\pm \mathrm{9.1}$ in the observations and the HadCM3 model respectively.

A3 SMB

Consistent with previous studies, the SMB is slightly too low in the Antarctic interior in HadCM3 (Turner et al., 2006), suggesting that the warm bias in these regions does not affect the modelled SMB. The largest model SMB errors (dry and wet biases) occur near the coasts (Fig. A2b). The dry biases may be due to the coarse HadCM3 grid, altering a realistic orography and representation of the ascending air masses that provide precipitation to these coastal regions. The coarse model grid biases can be seen in Fig. A2c as a step representation of the black points compared to an expected linear regression. Turner et al. (2006) also attribute the wet coastal biases to an overly intense mean sea level pressure field gradient: stronger air flows than in the observations produce excess precipitation on the western side of the Antarctic Peninsula. These aspects reduce the linear regression correlation (correlation coefficient of 0.70). The Antarctic mean climatological SMB difference between the observations and HadCM3 is −29.7 mm per month.

A4 δ¹⁸O

The distribution of the simulated δ¹⁸O over Antarctica is similar to observations (Antarctic means of $-\mathrm{36.2}\pm \mathrm{9.7}$ ‰ and $-\mathrm{37.4}\pm \mathrm{10.3}$ ‰ in the observations and the HadCM3 simulations respectively, minimum values of −61.3 ‰ and −57.9 ‰ in the observations and the HadCM3 simulations respectively, and maximum values of −3.2 ‰ and −7.7 ‰ in the observations and the HadCM3 simulations respectively). Excessively depleted values occur in the Antarctic interior (Fig. A3). These are associated with the warm bias. Overly enriched values are observed over the Peninsula and the Weddell Sea coast, consistently with the wet bias in these regions. Nevertheless, the HadCM3 historical simulations do capture the δ¹⁸O observations relatively well, as shown by the strong relationship between the outputs and the observations (correlation coefficient of 0.84 and slope of 0.90±0.02 ‰ ‰⁻¹).

Figure A1 Surface air temperature (SAT) evaluation: (a) map of the time-averaged HadCM3 SAT distribution over the Antarctic resulting from the ensemble mean for the historical period (in °C), (b) SAT difference between the time-averaged HadCM3 outputs from the ensemble mean for the historical period and the corresponding SAT observations (in °C), and (c) linear regression between the time-averaged HadCM3 outputs from the ensemble mean for the historical period and the corresponding SAT observations (black points). The red line is a 1:1 data–model slope.

Figure A2 Surface mass balance (SMB) evaluation: (a) map of the time-averaged HadCM3 precipitation minus evaporation (P−E) distribution over the Antarctic resulting from the ensemble mean for the historical period (in mm per month), (b) SMB difference between the time-averaged HadCM3 outputs from the ensemble mean for the historical period and the corresponding observations (in mm per month), and (c) linear regression between the time-averaged HadCM3 outputs from the ensemble mean for the historical period and the corresponding SMB observations (black points). The red line is a 1:1 data–model slope.

Figure A3 δ¹⁸O evaluation: (a) map of the time-averaged HadCM3 δ¹⁸O distribution over the Antarctic resulting from the ensemble mean for the historical period (in ‰), (b) SMB difference between the time-averaged HadCM3 outputs from the ensemble mean for the historical period and the corresponding observations (in ‰), and (c) linear regression between the time-averaged HadCM3 outputs from the ensemble mean for the historical period and the corresponding δ¹⁸O observations (black points). The red line is a 1:1 data–model slope.

Appendix B: Historical HadCM3/A2k δ¹⁸O comparison

Table B1 Regional surface air temperature trends (in °C per 100 years): A2k reconstructions over the last 100 years (Stenni et al., 2017), based on the ECHAM5-wiso model and scaled on the climate field reconstruction from Nicolas and Bromwich (2014) (“A2k lower bound” and “A2k upper bound”), and the HadCM3-simulated trends over the whole historical period (1851–2004). The relationships are significant.

Figure B1 Historical δ¹⁸O time series based on 5-year bins from ice core data and HadCM3. Time series of Antarctic surface-weighted averages of δ¹⁸O (in ‰) based on 5-year bin averages from ice core data (Stenni et al., 2017, solid black lines) and simulated by the HadCM3 model (ensemble mean, coloured lines) for each region of Antarctica as defined in the A2k project over the historical period 1851–2004. Dashed lines correspond to linear regressions.

Appendix C: Regional sea ice trends

Table C1 Regional sea ice area trends: slope plus or minus the standard error of the slope (in 10³ km² yr⁻¹). The correlation coefficients are given in parentheses. Only significant relationships are given. The sea ice regions are defined in terms of longitude as follows: the Indian sector is limited between 20 and 90° E. The Pacific sector is limited between 90 and 160° E. The Ross sector is limited between 160 and 230° E. The Bellingshausen–Amundsen sector is limited between 230 and 300° E. Finally, the Weddell sector is limited between 300 and 20° E.

Appendix D: Dependency of the time window on the HadCM3 δ¹⁸O trends

Table D1 Antarctic regional δ¹⁸O trends (Plateau, Indian Coast, Weddell Coast, Peninsula, WAIS, Victoria Land, and Dronning Maud Land, with region delimitations as defined in Stenni et al., 2017) simulated by the HadCM3 model using our historical mean ensemble using different time window lengths: the last 153, 100, 75, 60, 45, and 30 simulated years. The gradient and the standard error of the gradient are given for each linear relationship, associated with the correlation coefficient in brackets. Non-significant relationships are indicated by a dash.

Appendix E: Statistical description of the HadCM3 cold and warm ensembles

Table E1 Statistical description of the surface air temperature (SAT, in °C), precipitation (Precip, in mm per month), and δ¹⁸O (in ‰) for the cold and warm ensembles: minimum (Min), maximum (Max), and mean ± the standard deviation (mean ± SD).

Appendix F: δ¹⁸O–SAT gradients at ice core locations simulated by HadCM3

Table F1 δ¹⁸O–SAT gradients in ‰ °C⁻¹ at ice core locations. Gradients are accompanied by the standard error. Correlation coefficients are given in brackets. Finally, numbers in italics correspond to non-significant relationships (p>0.05).

Appendix G: Impact of the Southern Annular Mode

The Southern Annular Mode (SAM) is the leading mode of atmospheric variability in the Southern Hemisphere (Thompson and Wallace, 2000). In particular, it describes the position and the strength of the polar jet position, the southern westerly belt, and the associated storm tracks. A positive (negative) phase of the SAM is associated with an intensified (weakened) poleward (northward) shift of the polar jet. The SAM is thus the preferred studied mode to investigate the Southern Hemisphere teleconnection with lower latitudes.

Here, we used the definition of the SAM index following the approach of Gong and Wang (1999) as the difference between the normalised monthly zonal mean sea level pressure between 40 and 65° S. Here, we used the period 1961–1990 as a reference interval.

$$\begin{array}{}\text{(G1)}& \mathrm{SAM}={\displaystyle \frac{P_{\mathrm{40}}-{\mathit{\mu }}_{\mathrm{40}}}{{\mathit{\sigma }}_{\mathrm{40}}}}-{\displaystyle \frac{P_{\mathrm{65}}-{\mathit{\mu }}_{\mathrm{65}}}{{\mathit{\sigma }}_{\mathrm{65}}}},\end{array}$$

where P₄₀ and P₆₅ are the monthly mean sea level pressure at 40 and 65° S, μ₄₀ and μ₆₅ are the mean of the monthly mean sea level pressure at 40 and 65° S over the reference interval 1961–1990, and σ₄₀ and σ₆₅ are the standard deviations of the monthly mean sea level pressure at 40 and 65° S over the reference interval 1961–1990.

We computed the linear regressions between the calculated SAM and our climate variables (Fig. G1): (i) the surface air temperature (SAT), (ii) the precipitation (P), and (iii) the precipitation-weighted δ¹⁸O (δ¹⁸O). These linear regressions were computed over the whole historical simulated period and for the recent period 1950–2004 at the annual scale. Note that, as in the main text of the paper, we computed these linear regressions using the stack of the ensemble members, resulting in 918 points for the historical period (1851–2004) and 324 points for the period 1950–2004.

Figure G1 Correlation coefficients between the Southern Annular Mode index and the surface air temperature (SAT, a, d), the precipitation (P, b, e), and the precipitation-weighted δ¹⁸O (c, e) simulated by the HadCM3 model at the annual scale for the historical period (1851–2004, a–c) and the 1950–2004 period (d–f). Only significant relationships are shown (p-value < 0.05).

Within the frame of the CMIP5 project, the ability of HadCM3 to reproduce the SAM was evaluated (Zheng et al., 2013). As for all the CMIP5 models, HadCM3 overestimates the SAM index variability (Zheng et al., 2013; Zhang et al., 2022). Nevertheless, it reproduces the decadal variability of the SAM index and displays the best correlation coefficient between the modelled and observed detrended SAM index (Zheng et al., 2013).

Previous studies reported, based on observations, that the main Antarctic continent is globally colder and drier while the SAM is in a positive phase, as the stronger southern westerly wind belt reduces the exchanges with warmer air masses from midlatitude regions, with the exception of the Peninsula (Clem et al., 2016). These effects are reproduced in our HadCM3 simulations, as shown by the correlation coefficient values between the SAM and the SAT that are positive over the northern Antarctic Peninsula but negative over the rest of the continent, especially in coastal areas (Fig. G1a and d).

Similarly, it was shown that there is less southward moisture advection towards the Antarctic interior in a positive phase of the SAM, reducing precipitations. In our simulations (Fig. G1b, e), this effect is enhanced over the Antarctic Plateau, Victoria Land, and Marie Byrd Land. Conversely, the Antarctic Peninsula receives more precipitation. However, the discrepancy in the HadCM3 orography disables the “shadow effect”, decreasing precipitation on the eastern part of the Peninsula due to the presence of mountains (Fogt and Marshall, 2020).

The link between water stable isotopes and the SAM is less settled. A couple of publications displayed a correlation between the water stable isotope content in ice cores and the SAM index, but no systematic method allowed an established link. For instance, Servettaz et al. (2023a) suggest some impacts of the SAM on the isotopic content of the Aurora Basin North ice core over the last millennium, although not on the whole length of the core. Also, Vega et al. (2016) suggest that, over the Fimbul Ice Shelf, the absence of correspondence between water stable isotopes and SAT might be explained by changes in atmospheric circulation, supported by a high correlation between d-excess measured in the KM and BI ice cores and the SAM index. Kino et al. (2021) showed the contribution of SAM over precipitation-weighted δ¹⁸O at the daily scale simulated by the MIROC5-iso model nudged toward the JRA-25 reanalyses over the period 1981–2010 at Dome Fuji. However, they warn that it does not prevail in all Antarctic locations of the Antarctic Plateau. For instance, Dome C is less sensitive to SAM compared to possible other teleconnection modes (Dreossi et al., 2024). In our simulations, the correlation coefficients between the SAM and precipitation-weighted δ¹⁸O are significant and negative over the whole continent (Fig. G1c and f) but remain weak, with a mean of $-\mathrm{0.26}\pm \mathrm{0.11}$ over the historical period and $-\mathrm{0.27}\pm \mathrm{0.12}$ for the period 1950–2004. Thus, from our simulations, we cannot establish a robust link between the SAM and the Antarctic precipitation-weighted δ¹⁸O.

However, studying the impact of the atmospheric circulation change on Antarctic precipitation-weighted δ¹⁸O should not boil down to the link with the SAM. For instance, only some El Niño–Southern Oscillation (ENSO)/SAM combinations (El Niño/negative SAM and La Niña/positive SAM) contribute to strengthening the Amundsen Sea Low (e.g. Wilson et al., 2016), as observed through the analysis of the Roosevelt Island Climate Evolution (RICE) δ¹⁸O (Emanuelsson et al., 2023). SAM-induced processes impacting Antarctic precipitation-weighted δ¹⁸O are also not trivial: SAM changes SAT and precipitation regimes but also the sea ice in a more complex manner (Fogt and Marshall, 2020). Other modes affect the Antarctic atmospheric circulation and might explain the δ¹⁸O changes, as for the Indian Ocean Dipole in phase with El Niño through the production of atmospheric rivers (Shields et al., 2022).