Water stable isotope spatio-temporal variability in Antarctica in 1960–2013: observations and simulations from the ECHAM5-wiso atmospheric general circulation model

Polar ice core water isotope records are commonly used to infer past changes in Antarctic temperature, motivating an improved understanding and quantification of the temporal relationship between δ¹⁸O and temperature. This can be achieved using simulations performed by atmospheric general circulation models equipped with water stable isotopes. Here, we evaluate the skills of the high-resolution water-isotope-enabled atmospheric general circulation model ECHAM5-wiso (the European Centre Hamburg Model) nudged to European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis using simulations covering the period 1960–2013 over the Antarctic continent.

We compare model outputs with field data, first with a focus on regional climate variables and second on water stable isotopes, using our updated dataset of water stable isotope measurements from precipitation, snow, and firn–ice core samples. ECHAM5-wiso simulates a large increase in temperature from 1978 to 1979, possibly caused by a discontinuity in the European Reanalyses (ERA) linked to the assimilation of remote sensing data starting in 1979.

Although some model–data mismatches are observed, the (precipitation minus evaporation) outputs are found to be realistic products for surface mass balance. A warm model bias over central East Antarctica and a cold model bias over coastal regions explain first-order δ¹⁸O model biases by too-strong isotopic depletion on coastal areas and underestimated depletion inland. At the second order, despite these biases, ECHAM5-wiso correctly captures the observed spatial patterns of deuterium excess. The results of model–data comparisons for the inter-annual δ¹⁸O standard deviation differ when using precipitation or ice core data. Further studies should explore the importance of deposition and post-deposition processes affecting ice core signals and not resolved in the model.

These results build trust in the use of ECHAM5-wiso outputs to investigate the spatial, seasonal, and inter-annual δ¹⁸O–temperature relationships. We thus make the first Antarctica-wide synthesis of prior results. First, we show that local spatial or seasonal slopes are not a correct surrogate for inter-annual temporal slopes, leading to the conclusion that the same isotope–temperature slope cannot be applied for the climatic interpretation of Antarctic ice core for all timescales. Finally, we explore the phasing between the seasonal cycles of deuterium excess and δ¹⁸O as a source of information on changes in moisture sources affecting the δ¹⁸O–temperature relationship. The few available records and ECHAM5-wiso show different phase relationships in coastal, intermediate, and central regions.

This work evaluates the use of the ECHAM5-wiso model as a tool for the investigation of water stable isotopes in Antarctic precipitation and calls for extended studies to improve our understanding of such proxies.

1 Introduction

The Antarctic climate has been monitored from sparse weather stations providing instrumental records starting at best in 1957 (Nicolas and Bromwich, 2014). Water stable isotopes in Antarctic ice cores are key to expanding the documentation of spatio-temporal changes in polar climate and the hydrologic cycle (Jouzel et al., 1997) for the recent past (PAGES 2k Consortium, 2013; Stenni et al., 2017) as well as for glacial–interglacial variations (Jouzel et al., 2007; Schoenemann et al., 2014). Water stable isotopes measured along ice cores were initially used to infer Antarctic past temperatures using the spatial isotope–temperature slope (Lorius et al., 1969). The focus on inter-annual variations is motivated by the goal of quantifying temperature changes at the Earth's surface, including Antarctica, during the last millennia, to place current changes in the perspective of recent natural climate variability (Jones et al., 2016), to understand the drivers of this variability, and to test the ability of climate models to correctly represent it. This timescale is relevant for the response of the Antarctic climate to e.g. volcanic forcing and for the Antarctic climate fingerprint of large-scale modes of variability such as ENSO and the Southern Annular Mode (Smith and Stearns, 1993; Turner, 2004; Stammerjohn et al., 2008; Schroeter et al., 2017). The various climate signals potentially recorded in precipitation isotopic composition are, however, difficult to disentangle.

First, the original signal from precipitation may be altered due to deposition and post-deposition processes (e.g. Sokratov and Golubev, 2009; Jones et al., 2017; Münch et al., 2017; Laepple et al., 2018). Wind erosion and sublimation during or after precipitation have long been known to affect ice core records (Eisen et al., 2008; Grazioli et al., 2017). Other processes such as melt and diffusion can also alter the preservation of isotopic signals in firn and ice and cause smoothing of the initial snowfall signals (Johnsen, 1977; Whillans and Grootes, 1985; Johnsen et al., 2000; Jones et al., 2017). So far, the mechanisms of such post-deposition processes on the alteration of the initial precipitation signals are not fully understood and quantified (Touzeau et al., 2017). Second, the Antarctic snowfall isotopic composition may be affected by the origin of moisture and the associated evaporation conditions, changes in the relationship between condensation and surface temperature, or changes in the intermittency of precipitation (e.g. Sime et al., 2009; Hoshina et al., 2014; Touzeau et al., 2016). Although the surface snow isotopic composition signal has classically been interpreted as a precipitation-weighted deposition signal (Krinner and Werner, 2003), recent studies evidenced isotopic exchanges between the Antarctic snow surface and the atmosphere associated with snow metamorphism occurring at diurnal and sub-annual scales (Steen-Larsen et al., 2014; Casado et al., 2016; Ritter et al., 2016; Touzeau et al., 2016).

Second, the climatic interpretation of water stable isotopes in Antarctic ice cores is still challenging. Quantitative approaches have relied on empirical relationships and the use of theoretical and atmospheric models including water stable isotopes. Pioneer studies evidenced a close linear relationship between the spatial distribution of water stable isotopes and local temperature (e.g. Lorius and Merlivat, 1975) and explained this feature as the result of the distillation along air mass trajectories. Thereupon, local temperature (i.e. at a specific site) was reconstructed using δ¹⁸O measurements and based on the slope of the aforementioned spatial empirical relationship as a surrogate for relationships at annual to multi-annual scales. However, recent data syntheses have shown that other effects had to be taken into account (e.g. Masson-Delmotte et al., 2008). It was found that Antarctic snowfall isotopic composition is also linked to the initial vapour isotopic composition (Stenni et al., 2016), atmospheric transport pathways (Schlosser et al., 2008; Dittmann et al., 2016), Antarctic sea ice extent (Bromwich and Weaver, 1983; Noone and Simmonds, 2004; Holloway et al., 2016), and local condensation temperature, which is itself related to surface temperature through complex boundary layer processes (Krinner et al., 2007). Evaporation conditions, transport, and boundary layer processes may vary through time from seasonal (Fernandoy et al., 2018) to annual or multi-annual scale, thereby potentially distorting the quantitative relationship between snow isotopic composition and local surface air temperature estimated empirically for present-day conditions (Jouzel et al., 1997).

Model studies have been key to quantitatively exploring the spatio-temporal aspects of the relationships between precipitation isotopic composition and temperature (Jouzel et al., 2000). Mixed-cloud isotopic models have been used to propose a coherent interpretation of δ¹⁸O and δD data in terms of changes in site and source temperatures (Uemura et al., 2012) or to simulate isotopic variations along individual atmospheric trajectories (Dittmann et al., 2016). However, such theoretical distillation models rely on the closure assumption at the ocean surface to calculate the initial evaporation isotopic composition and do not account for atmospheric dynamics and mixing of air masses (Jouzel and Koster, 1996; Delmotte et al., 2000). Atmospheric general circulation models equipped with water stable isotopes offer a physically coherent, three-dimensional framework to investigate the weather and climate drivers of Antarctic precipitation isotopic composition (Jouzel et al., 2000). They play a key role in assessing how different boundary conditions (e.g. changes in orbital forcing, changes in atmospheric greenhouse gas concentration) affect the simulated relationships between precipitation isotopic composition and climate variables. Most of these simulations support the idea that the present-day isotope–temperature spatial relationship is a good approximation for the relationships between glacial conditions and today (Delaygue et al., 2000; Werner et al., 2018), with one exception (Lee et al., 2008). One study used climate projections in response to increased atmospheric CO₂ concentration to explore isotope–temperature relationships in a world warmer than today and suggested a changing temporal isotope–temperature relationship due to changing covariance between temperature and precipitation (Sime et al., 2009). Several observational and modelling studies have also evidenced different isotope–temperature relationships between the spatial relationship and those calculated at the seasonal (Morgan and van Ommen, 1997) or inter-annual scale (Schmidt et al., 2007).

Our study is motivated by the need for a synthesis over all of Antarctica using a proper interpretation of processes that affect water stable isotopes on the appropriate spatial and temporal scales. It aims to address the following questions: (i) what is the performance of a state-of-the-art atmospheric general circulation model with respect to existing Antarctic observations of spatio-temporal variations in temperature, surface mass balance, precipitation, and snow isotopic composition for present day? (ii) What can we learn from such a model for the regional relationships between isotopic composition from the precipitation and temperature at the inter-annual scale for the recent past and considering all of Antarctica?

Figure 1 Spatial distribution of Antarctica in seven regions: East Antarctic Plateau, coastal Indian, Weddell sea, West Antarctic Ice Sheet, Victoria Land, and Dronning Maud Land. The location of the selected READER surface stations: Neumayer, Mawson, Vostok, Dome C, Casey, Dumont d'Urville (noted as “DDU”), McMurdo, Byrd, Palmer, and Esperanza.

For this purpose, we focus on the high-resolution atmospheric general circulation model equipped with water stable isotopes, ECHAM5-wiso (the European Centre Hamburg model), which demonstrated remarkable skills for Antarctica (Werner et al., 2011). We explore a simulation performed for the period 1960–2013 in which the atmospheric model is nudged to the European Reanalyses (ERA) ERA-40 and ERA-Interim (Uppala et al., 2005), ensuring that the day-to-day simulated variations are coherent with the observed day-to-day variations in synoptic weather and atmospheric circulation (see Butzin et al., 2014, for more explanation). This framework is crucial to performing comparisons between simulations and observations for temporal variations. Second, we compile a database of precipitation, snow, and firn–ice isotopic composition using data from precipitation sampling and ice core records and considering δ¹⁸O and deuterium excess (hereafter, d). These methods are described in Sect. 2. We then compare the model outputs with the available datasets (Sect. 3). After evaluating the near-surface temperature and the surface mass balance (hereafter SMB; Sect. 3.1), we focus on the water stable isotopes (Sect. 3.2). We emphasise spatial patterns, the magnitude of inter-annual variability (Sect. 3.2.1 and 3.2.4), and the pattern and the amplitude of seasonal variations (Sect. 3.2.2 and 3.2.4). We explore the simulated and estimated isotope–temperature relationships (Sect. 3.2.3) and the relationships between d and δ¹⁸O (Sect. 3.2.3). Highlighting the strengths and limitations of the model (Sect. 3.3), we use the simulation framework to explore the δ¹⁸O–temperature relationship (Sect. 4.1) and the phase lag between seasonal variations in d and δ¹⁸O (Sect. 4.2). Finally, we focus on the implications of our results for the climatic interpretation of water stable isotope records for seven Antarctic regions (central plateau, coastal Indian, Weddell Sea coast, West Antarctic Ice Sheet, Victoria Land, and Dronning Maud Land). The Antarctica2k group (Stenni et al., 2017) indeed identified these seven Antarctic regions, which are geographically and climatically consistent, to produce regional temperature reconstructions using ice core records. The results of our study thus contribute to the reconstruction of past Antarctic climate spanning the last 2000 years (the Antarctica2k initiative) of the Past Global Changes (PAGES) PAGES2K project (PAGES 2k Consortium, 2013) by providing quantitative calibrations of the regional temperature reconstructions using ice core water stable isotope records.

2 Material and methods

2.1 Observations and reanalysis products

2.1.1 Temperature and surface mass balance instrumental records

Station temperature records have been extracted from the READER database (https://legacy.bas.ac.uk/met/READER, last access: August 2017; Turner et al., 2004). We have selected surface station data following two conditions: to cover the seven Antarctic regions aforementioned (see Sect. 1 and Fig. 1) with at least one station for each and to cover the period 1960–2013. As a result, we have selected Neumayer, Mawson, Vostok, Casey, Dumont d'Urville (hereafter DDU), McMurdo, Palmer, and Esperanza station surface data. Due to the short duration of surface station records for the 90–180^∘ W sector, we have added data from the automatic weather station (hereafter, AWS) of Dome C, but we have used it with caution as these records are associated with a warm bias in thermistor measurements due to solar radiation when the wind speed is low (Genthon et al., 2011). Finally, we extracted the reconstruction of temperature for Byrd station by Bromwich et al. (2013) based on AWS data and infilled with observational reanalysis data. No record meets our criteria for the Weddell Sea coast region (Fig. 1).

SMB data have been extracted from the quality-controlled GLACIOCLIM-SAMBA (GC) database (Favier et al., 2013). We have selected data spanning the twentieth century, corresponding to 3242 punctual values, which have then been clustered within the corresponding ECHAM5-wiso grid cells for the calculation of gridded annual average values. As described by Favier et al. (2013), the spatial coverage of SMB field data is particularly poor in the Antarctic Peninsula, in West Antarctica, and along the margins of the ice sheet. As a result, SMB is not correctly sampled at elevations between 200 and 1000 m a.s.l., where accumulation rates are the highest. In central Antarctica, areas characterised by wind glaze and megadunes are also insufficiently documented.

2.1.2 ERA reanalyses

The ECHAM5-wiso model run for this study is nudged to ERA-40 (Uppala et al., 2005) and ERA-Interim (Dee et al., 2011) global atmospheric reanalyses produced by the European Centre for Medium-Range Weather Forecasts (ECMWF). ERA-40 covers the period 1957–2002 at a daily resolution, with a spatial resolution of 125 km × 125 km. ERA-Interim covers the period 1979 to present at a 6-hourly resolution, with a spatial resolution of 0.75^∘ × 0.75^∘.

For comparison with instrumental records and ECHAM5-wiso outputs, we have extracted 2 m temperature outputs (hereafter 2 m-T) over the periods 1960–1978 and 1979–2013 for ERA-40 and ERA-Interim, respectively, at grid cells closest to the stations where meteorological measurements have been selected (see previous section). We have then calculated annual averages.

2.1.3 A database of Antarctic water stable isotopic composition from precipitation, surface snow, and firn–ice core records

This database consists of water stable isotope measurements performed on different types of samples (precipitation, surface snow, or shallow ice cores) and at different time resolutions (sub-annual, annual, or multi-annual average values; see Table S1 in the Supplement). Sample data consist of δ¹⁸O and/or δD, providing d, if both δ¹⁸O and δD have been measured. Altogether, we have gathered data from the following.

1. A total of 101 high-resolution ice core records, including 79 annually resolved records and 18 records with sub-annual resolution (including 5 records with both δ¹⁸O and δD data). These data have been extracted from the Antarctica2k data synthesis (Stenni et al., 2017) with a filter for records spanning the interval 1979–2013, thus restricting the original 122 ice cores to a resulting 101 ice core data. Primary data sources, geographical coordinates, and covered periods are reported in Table S1 in the Supplement.

2. Average surface snow isotopic composition data compiled by Masson-Delmotte et al. (2008; available at http://www.lsce.ipsl.fr/Phocea/Pisp/index.php?nom=valerie.masson, last access: August 2017) expanded with datasets from Fernandoy et al. (2012); in this case, the averaging period is based on different time periods, with potential non-continuous records (see Table S1 in the Supplement).

3. Precipitation records extracted from the International Atomic Energy Agency/Global Network of Isotopes in Precipitation (IAEA/GNIP) network (IAEA/WMO, 2016), with monthly records available for four Antarctic stations, complemented by daily records for four Antarctic stations from individual studies. Precipitation records from Vostok are available but are excluded from our analysis due to an insufficient number of measurements (29). See orange part of Table S1 in the Supplement.

Each of the 1205 locations have given an individual index number. Data have been processed to calculate time-averaged values (available at 1089 locations for δ¹⁸O values, 879 locations for δD, and 770 locations for d). The ice core records with sub-annual resolution were averaged at annual resolution over the period 1979–2013, resulting in 88 ice core records for δ¹⁸O and only 5 for d. Most precipitation records are not continuous and do not cover a full year, preventing the calculation of annual mean values. We have also used sub-annual records from 22 highly resolved ice cores (including 18 records giving access to δ¹⁸O and 5 records giving access to d) and precipitation sampling from eight stations to characterise the seasonal amplitude. For ice core records, we have only calculated the yearly amplitude from available measurements, as chronologies cannot be established at monthly scales. Note that this database is publicly available on the PANGAEA data archive (https://www.pangaea.de/?t=Cryosphere, last access: June 2018).

2.2 ECHAM5-wiso model and simulation

The atmospheric general circulation model (AGCM) ECHAM5-wiso (Roeckner et al., 2003; Werner et al., 2011) captures the global pattern of precipitation and vapour isotopic composition, including the spatial distribution of annual mean precipitation isotopic composition over Antarctica (Masson-Delmotte et al., 2008). Several studies using ECHAM5-wiso have been dedicated to model–data comparisons for temporal variations in other regions (e.g. Siberia, Greenland; Butzin et al., 2014; Steen-Larsen et al., 2016).

The ECHAM5-wiso outputs analysed in this study consist of daily values simulated over the period 1960–2013. ECHAM5-wiso was nudged to atmospheric reanalyses from ERA-40 (Uppala et al., 2005) and ERA-Interim (Dee et al., 2011), which are shown to have good skills for Antarctic precipitation (Wang et al., 2016), surface pressure fields, and vertical profiles of winds and temperatures. The ocean surface boundary conditions (sea ice included) are also prescribed based on ERA-40 and ERA-Interim data. Isotope values of ocean surface isotopic composition are based on a compilation of observational data (Schmidt et al., 2007). The simulation was performed at a T106 resolution (which corresponds to a mean horizontal grid resolution of approx. 1.1^∘ × 1.1^∘) with 31 vertical model levels.

2.3 Methods for model–data comparisons

In the model, we have extracted specific daily variables for comparison with available data and then averaged them. We have extracted daily 2 m temperature outputs (hereafter 2 m-T) for comparison with surface air instrumental records, daily (precipitation minus evaporation) outputs (hereafter P-E) for comparison with SMB data, and daily precipitation isotopic composition outputs for comparison with measurements of isotopic composition data in the precipitation. For ice core data, we averaged daily precipitation isotopic composition weighted by the daily amount of precipitation.

For each specific site, we selected the model grid cell including the coordinates of the site. When comparing model outputs with the database of surface data (time-averaged SMB and isotopic composition), available data have been averaged within each model grid cell.

Time selection was dependent on the variables. The 2 m-T outputs have been compared with temperature records for the period 1960–2013 based on annual averages and selecting the same years as in the data (see Sect. 3.1.1). The comparison with other datasets (SMB, snow, and water stable isotopes from firn–ice cores) is restricted to the period 1979–2013 due to concerns about the skills of the reanalyses used for the nudging prior to 1979 in Antarctica (see next section). Daily (P-E) outputs were all extracted over the whole period 1979–2013 and averaged (see Sect. 3.1.2). For comparison with the surface isotopic database (Sect. 3.2.1), daily precipitation isotopic composition was averaged by weighting by the daily amount of precipitation over the whole period 1979–2013. For the inter-annual variability (same section) or annual values (e.g. for d outputs, see Sect. 4), daily precipitation isotopic composition was averaged by weighting by the daily amount of precipitation for each year of the period 1979–2013. For sub-annual isotopic composition, we used precipitation isotopic compositions (amplitude and mean seasonal cycle) and highly resolved ice cores (amplitude only). Precipitation isotopic composition data consist of a very small number of measurements, sometimes taken before 1979 (e.g. observations from DDU consist of 19 measurements during 1973), and thus model precipitation isotopic composition outputs were extracted at the very exact sampling date. Then, monthly averages were performed and mean seasonal cycles were calculated. The resulting mean seasonal cycles of precipitation isotopic composition were obtained the same way in both precipitation data and the model. For comparison with the mean seasonal amplitude of the highly resolved ice cores, the mean seasonal amplitude was calculated from the mean seasonal cycle based on the monthly averages (weighted by the precipitation amount) over the period covered by the ice core record.

Finally, for the spatial linear relationships, the calculations reported for each grid cell are based on the relationship calculated by including the 24 grid cells (±2 latitude steps, ±2 longitude points) surrounding the considered grid cell.

Our comparisons are mainly based on linear regressions. Note that through the paper, we consider a linear relationship to be significant for a p value < 0.05.

3 Model skills

Figure 2 Surface air temperature (in ^∘C) from station instrumental records (points and dashed lines) and simulated by the ECHAM5-wiso model (solid lines) over the period 1960–2013 for (a) the plateau, (b) coastal East Antarctic Ice Sheet, and (c) the West Antarctic Ice Sheet. Note that the plots were organised by regions to make it more readable: inland (a), coastal (b), and West Antarctic Ice Sheet plus peninsula (c).

Figure 3 The 2 m temperature outputs (in ^∘C) from ERA-40 (light green), ERA-Interim (dark green), and ECHAM5-wiso outputs over the periods 1960–1979 (light purple) and 1979–2013 (dark purple) at the locations of Neumayer, Byrd, Palmer, Vostok, Dome C, McMurdo, Casey, Dumont d'Urville (written as DDU), Mawson, and Esperanza stations. Horizontal black lines correspond to the mean data. Vertical black lines correspond to inter-annual standard deviations: dotted lines are associated with data, while solid lines are associated with model outputs (ERA or ECHAM).

In this section, we assess ECHAM5-wiso skills with the perspective of using the model outputs for the interpretation of water stable isotope data. In polar regions, isotopic distillation is driven by fractionation occurring during condensation, which is itself controlled by condensation temperature (Dansgaard, 1964). We thus first compare ECHAM5-wiso outputs with regional climate records, as this comparison may explain potential isotopic biases. This includes a comparison with reanalyses in order to explore the role of nudging in model–data mismatches. We then compare ECHAM5-wiso outputs with our isotopic database.

3.1 Temperature and surface mass balance

3.1.1 Comparison with instrumental temperatures records and ERA outputs

Table 1 Differences between observed (READER) and simulated (ECHAM5-wiso) annual surface air temperature: observed average (noted as “observed μ”, in ^∘C), average difference (noted as “μ differences”, in  ^∘C), standard deviation from observations (noted as “observed σ”, in ^∘C), and standard deviation from the model (noted as “simulated σ”, in ^∘C) for the period 1979–2013.

We compare time series of instrumental temperature records (filled circles and dashed lines, Fig. 2) with model outputs (solid lines, Fig. 2) from 1960 to 2013. This comparison first highlights local offsets between observed and simulated mean values at each site, without a systematic overall warm or cold bias. Table 1 reports the statistical analysis of annual differences between observations and simulations (observed mean, mean difference between the data and the model outputs, observed versus simulated standard deviation). ECHAM5-wiso has a cold bias for 7 out of 10 stations. While this bias is less than 2 ^∘C for Dronning Maud Land (Mawson and Neumayer) and over the peninsula (Palmer and Esperanza), it reaches 7 ^∘C for the coastal Indian region (Casey and Dumont d'Urville) and is very strong over the Victoria Land region (McMurdo), reaching 15 ^∘C. This cold bias may be due to the model resolution and the location of coastal stations in the ice-free region, where the small-scale topographic features are not accounted for at the model resolution. In contrast, ECHAM5-wiso has a warm bias for all the stations located inland (Vostok, Dome C, and Byrd). Werner et al. (2011) also reported this warm bias for the central Antarctic Plateau and suggested that it could be linked to problems in correctly simulating the polar atmospheric boundary layer. Our comparison also shows that the simulated inter-annual temperature variability is larger than observed for seven out of nine sites and is particularly overestimated for locations such as DDU, McMurdo, and Palmer, where the cold bias is large.

Figure 2 depicts a sharp simulated increase in temperature from 1978 to 1979 for all stations, except for the peninsula region (Esperanza and Palmer). Such a feature is not displayed in instrumental records, with one exception at McMurdo (Fig. 2). As a result, the model–data correlation coefficient for McMurdo is higher over 1960–2013 than over 1979–2013 (Table 2), possibly because it is dominated by the sharp increase just prior to 1979. For all other stations, the correlation coefficient is significantly higher in 1979–2013 than in 1960–2013. In order to assess whether ECHAM5-wiso reproduces the temperature bias displayed by ERA-40 (Bromwich et al., 2007), we compare outputs from ERA-40 and ERA-Interim (green bars, Fig. 3) with ECHAM5-wiso outputs (purple bars, Fig. 3) nudged by these reanalyses (i.e. over 1960–1978 and 1979–2013, respectively) and with the station temperature data (horizontal black lines, Fig. 3).

Table 2 Linear relationship between surface temperatures (in ^∘C) from station instrumental records and ECHAM5-wiso outputs (in ^∘C) over the period 1960–2013 and 1979–2013: the slope (in ^∘C ^∘C⁻¹), the correlation coefficient (noted as “r”), and the p value. Data are not reported for 1960–2013 for stations for which records only cover the second period (1979–2013). Numbers in brackets correspond to standard errors.

All datasets reveal a cold bias simulated by both the reanalyses and ECHAM5-wiso at all stations but Byrd and Vostok over the two periods (only over 1960–1978 for Neumayer and Esperanza), but this bias is larger over the period 1960–1978 compared to the period 1979–2013. This finding supports our earlier suggestion for Dumont d'Urville (Goursaud et al., 2017) that the 1978–1979 shift simulated by ECHAM5-wiso arises from the nudging to ERA-40 reanalyses. We note that mean values and the amplitude of inter-annual variations are different for ECHAM5-wiso and ERA (not shown), as expected from different model physics despite the nudging technique. This finding has led us to restrict, as far as possible, the subsequent analysis of the ECHAM5-wiso outputs to the period 1979–2013.

Figure 4 Comparison of the GLACIOCLIM (hereafter and noted in the plots as “GC”) SMB database averaged within the ECHAM5-wiso grid cells and the SMB (i.e. precipitation–evaporation) simulated by the model, with the spatial distribution of the accumulation (a) as simulated by the model (in cm w.e. yr⁻¹), (b) the ratio of the ECHAM5-wiso annual accumulation (precipitation minus evaporation) to the GC averaged SMB (no unit), and GC-averaged SMB values against SMB values simulated by the model (blue dots) associated with the corresponding linear relationships (red solid line); displayed at the logarithm scale for elevation ranges of 0–2200 m a.s.l. (with the upper limit excluded) (c) and 2200–4000 m a.s.l. (d).

Figure 5 Maps displaying model–data comparisons for δ¹⁸O time-averaged values (a) and inter-annual standard deviations (b). Backgrounds correspond to ECHAM5-wiso simulations over the period 1979–2013, while signs correspond to the model–data comparison. For the time-averaged values, the comparison consists of calculating the model–data differences. Red “+” symbols indicate a positive model–data difference, while blue “−” symbols correspond to a negative model–data difference. For the inter-annual standard deviations, the comparison consists of calculating the ratio of the simulated value to the corresponding grid cell data. Red “+” symbols indicate a ratio higher than 1, while blue “−” symbols correspond to a model ∕ data ratio lower than 1.

Figure 6 Maps displaying model–data comparisons for d time-averaged (in ‰, a) values and inter-annual standard deviations (in ‰, b). Backgrounds correspond to ECHAM5-wiso simulations over the period 1979–2013, while signs correspond to the model–data comparison. For the time-averaged values, the comparison consists of calculating the model–data differences. Red “+” symbols indicate a positive model–data difference, while blue “−” symbols correspond to a negative model–data difference. For the inter-annual standard deviation, the comparison consists of calculating the ratio of the simulated value to the corresponding grid point data. Red “+” symbols indicate a ratio higher than 1, while blue “−” symbols correspond to a model ∕ data ratio lower than 1.

For this period marked by small temperature variations, we note that the correlation coefficient between data and model outputs (Table 2) is very small for McMurdo (r = 0.2) and rather small for Vostok (r = 0.6), questioning the ability of our simulation to resolve the drivers of inter-annual temperature variability at these locations. We observe that the model reproduces the amplitude of inter-annual variations, with a tendency to underestimate the variations as shown by model–data slopes from 0.6 to 1 ^∘C per ^∘C. As a result, ECHAM5-wiso underestimates the magnitude of inter-annual temperature variability for these central regions of the West and East Antarctic Ice Sheet. It will therefore be important to test whether similar caveats arise for water isotopes.

3.1.2 Comparison with GLACIOCLIM database accumulation

For each grid cell in which at least one stake record is available, we have calculated the ratio of the P-E values (which we use as a surrogate for accumulation) simulated by ECHAM5-wiso to the averaged SMB estimate for that grid region based on stake measurements (Fig. 4a). Due to the limited number of grid cells containing SMB data points from 1979 to 2013 (100 cells) located almost only on the East Antarctic Ice Sheet, we have decided to use the dataset covering the entire twentieth century (521 cells) spread over the continent.

The spatial distribution of SMB is well captured by ECHAM5-wiso, with decreasing SMB values from the coast to the interior plateau (Fig. 4a). However, the model quantitatively shows some discrepancies when compared with the GC database. The area-weighted (by the model grid cells) mean GC SMB is 141.3 mm w.e. yr⁻¹, while the simulated area-weighted mean P-E over the same model grid is 126.6 mm w.e. yr⁻¹. This underestimation covers 69.7 % of the compared areas. The 30.3 % remaining areas associated with an overestimation of the model are located in sparse regions like in the north of the plateau and over coastal areas (Fig. 4b). Note that the low P-E rates over the plateau (75 mm w.e. yr⁻¹, see Fig. 4a) counterbalance the local overestimation at the coast, supporting the ability of ECHAM5-wiso to resolve the integrated surface mass balance for the Antarctic ice sheet. Figure 4c and d confirm the global underestimation by the model, with slopes of simulated P-E against GC SMB lower than 1. This aspect is emphasised for elevations higher than 2200 m a.s.l. (r = 0.74 and rmse = 122.8 mm w.e. yr⁻¹ for elevation lower than 2200 m a.s.l., and r = 0.83 and rmse = 55 mm w.e. yr⁻¹ for elevations higher than 2200 m a.s.l., with “r” the correlation coefficient, and “RMSE” the root mean square error). The correlation coefficient (considering all elevations) is 0.79, reflecting the non-homogenous bias over the whole continent. This can be due first to a failure in the representativity of SMB spatial variability when averaging GC data within ECHAM5-wiso grid cells due to a too-small number of point measurements. Second, the model grid resolution may be too coarse to reproduce coastal topography and thus the associated amounts of precipitation. Finally, several key processes such as blowing snow erosion and deposition are not taken into account in the model. For instance, the lowest value from the GC database is −164 mm w.e. yr⁻¹, measured at the Bahia del Diablo glacier, a small glacier covering important elevation ranges in a narrow spatial scale between the front and the summit. It was the only one within the corresponding model grid cell, so the resulting GC value within this grid cell could not be representative of the model scale, and the simulated P-E value is not representative of this small glacier-wide value.

When considering the whole Antarctic grounded ice sheet, the area-weighted P-E simulated by the model amounts to 164.4 mm w.e. yr⁻¹. This value falls within the highest values of the 11 simulations displayed by Monaghan et al. (2006), varying from 84 to 188 mm w.e. yr⁻¹. However, the high range of values between the different simulations illustrates the uncertainties related to the SMB model, mainly due to model resolution, which is crucial to reproducing the impact of topography on precipitations and to non-resolved physical processes (e.g. drifting snow transport, including the erosion, deposition, and sublimation of drifting snow particles, and clouds microphysics; Favier et al., 2018). Moreover, this simulated value is very close to the best estimations of Antarctic grounded ice sheet SMB, which range between 143.4 (Arthern et al., 2006) and 160.8 mm w.e. yr⁻¹ (Lenaerts et al., 2012). This simulated value is also very close to the one obtained by Agosta et al. (2013) for the LMDZ4 model over the period 1981–2000 (160 mm w.e. yr⁻¹), but slightly lower than with the SMHiL model forced by LMDZ4 (189 mm w.e. yr⁻¹).

Table 3 Comparison between measurements from precipitation samples and ECHAM5-wiso simulated precipitation isotopic composition for grid cells closest to the sampling locations over the same period as the data (at daily or monthly scale when the name of the station is associated with an asterisk). We report the mean value ± the standard deviation for δ¹⁸O (in ‰) and for temperature (^∘C).

To conclude, although the ECHAM5-wiso simulation presented in this study has a relatively coarse resolution (110 km × 110 km compared to 15 km × 15 km for the SMHiL model forced by LMDZ4) and does not resolve processes contributing in the SMB (e.g. drifting snow processes), the P-E outputs are realistic products when compared with SMB data.

3.2 Comparison with water stable isotope data

Limited by the availability of the data, we could only study model skills with respect to spatio-temporal patterns, including seasonal and inter-annual variations, and the simulated relationships between δ¹⁸O and temperature. We have also extended the model–data comparison to the second-order parameter, d.

3.2.1 δ¹⁸O time-averaged values and inter-annual variability

The model–data difference of the time-averaged values is positive for 88 % of all grid cells, suggesting a systematic underestimation of isotopic depletion by ECHAM5-wiso (Fig. 5a). The few areas for which ECHAM5-wiso overestimates the isotopic depletion are restricted to coastal regions. This pattern is coherent with the temperature anomalies: ECHAM5-wiso produces too-low isotopic values where ECHAM5-wiso has a cold bias, likely causing too-strong distillation towards coastal areas and too-high isotopic values inland, where the warm bias limits the distillation strength. The statistical distribution of model–data δ¹⁸O differences (not shown) shows a wide range but an interquartile range (50 % of all values) of 1.4 to 3.9 ‰, which is therefore within 1.3 ‰ of the median. We conclude that, beyond the systematic offset linked to climatic biases, ECHAM5-wiso correctly captures the spatial gradient (continental effect) of annually averaged δ¹⁸O data. These results also suggest that the spatial distribution of annual mean δ¹⁸O values from shallow ice cores is driven by transport and condensation processes well resolved by ECHAM5-wiso, probably with secondary effects of non-resolved processes such as snow drift, wind erosion, and snow metamorphism. The largest deviations are encountered in coastal regions, where the model resolution is too low to correctly resolve topography, advection, and boundary layer processes (e.g. small-scale storms, katabatic winds). Katabatic winds also have the potential to enhance ventilation-driven post-deposition processes (Waddington et al., 2002; Neumann and Waddington, 2004).

δ¹⁸O inter-annual standard deviation is underestimated by the model for 92 % of the 179 grid cells in which this comparison can be performed (Fig. 5b).The interquartile range of the ratio between the simulated and observed standard deviation varies from 0.4 to 0.6 (not shown), with an underestimation by a factor of 2 for about 50 % of the grid cells. No such underestimation of inter-annual standard deviation was identified for the simulated temperature.

We now focus on our model–data comparison of precipitation data. Both precipitation isotopic composition and temperature measurements are available for only eight locations and for short time periods (Table 3). These data evidence the altitude and continental effect with increased isotopic depletion from Vernadsky (averaged δ¹⁸O of −9.9 ‰) to Dome F (averaged δ¹⁸O of −61.3 ‰). For five out of the eight records, the isotopic depletion is stronger in ECHAM5-wiso than observed (Dome C included). The observations depict an enhanced inter-daily δ¹⁸O standard deviation for inland sites, from 3.1 at Vernadsky to 10.8‰ at Dome F. The simulated δ¹⁸O inter-daily standard deviation is 1.1 to 3.8 times larger than observed, ranging from 5.1 to 19.2 ‰. For the exact same time period corresponding to the short precipitation isotopic records, ECHAM5-wiso simulates colder than observed temperatures at all stations but Dome F and Dome C, i.e. over the plateau. This finding is consistent with results from ice core records reported previously and consistent with the isotopic systematic biases. From this limited precipitation dataset, there is no systematic relationship between model biases for temperature (mean value or standard deviation) or for δ¹⁸O in contrast with the outcomes of the model–data comparison using the whole dataset, including surface snow. At Dome C, ECHAM5-wiso underestimates the standard deviation of temperature, but strongly overestimates the standard deviation of δ¹⁸O.

As a conclusion, while δ¹⁸O time-averaged model–data biases are consistent with temperature biases using the whole dataset, no systematic relationships emerge between model biases for temperature and δ¹⁸O measured in precipitation.

3.2.2 δ¹⁸O seasonal amplitude

High-resolution δ¹⁸O data allow us to explore seasonal variations. This includes 18 ice core records with sub-annual resolution, four IAEA/GNIP monthly precipitation datasets, and four daily precipitation monitoring records.

In order to quantify post-deposition effects in ice cores, we calculated the ratio of the first three seasonal amplitudes by using the mean seasonal amplitude in sub-annual ice cores (See Table S2 in the Supplement). We find a mean ratio of 1.40 ± 0.47. We explored whether this ratio was related to annual accumulation rates (see Fig. S3 in the Supplement), without any straightforward conclusion. We also observe that five ice cores depict a ratio lower than 1, including one with a mean yearly accumulation of 15 cm w.e. yr⁻¹, a feature which may arise from inter-annual variability in the precipitation seasonal amplitude or in post-deposition processes. This empirical analysis shows that a loss of seasonal amplitude due to post-deposition processes is likely in most cases, with an average loss of the seasonal amplitude of approximately 70 % compared to the amplitude recorded in the upper part of the firn cores (first 3 years).

We have calculated the mean of the δ¹⁸O annual amplitude (i.e. maximum–minimum values within each year) in ice core records (triangles in Fig. 7a) and the mean seasonal amplitude of precipitation time series (circles in Fig. 7a) for comparison with ECHAM5-wiso outputs (Fig. 7, Table 4). Unfortunately, a too-small number of measurements (19 daily measurements) were monitored at DDU, preventing the representation of the full seasonal cycle. The data depict the largest seasonal amplitude in the central Antarctic Plateau, reaching up to 25.9 ‰ at Dome F. ECHAM5-wiso underestimates the seasonal amplitude (by 14 to 69 %) when compared to precipitation data, but overestimates the seasonal amplitude when compared to ice core data (from 11 to 71 %). The overestimation when comparing with ice core data is consistent with the attenuation of signal by post-deposition effects (as previously mentioned) rather than a model bias.

Figure 7 Average seasonal amplitude of precipitation δ¹⁸O (a) and d (b) (in ‰) simulated by ECHAM5-wiso (colour shading) over the period 1979–2013 and calculated from precipitation data (circles) and ice core records (triangles) over their respective available periods.

Table 4 δ¹⁸O mean seasonal amplitude (in ‰) calculated for precipitation and sub-annual ice core data, as well as simulated by ECHAM5-wiso for the same time period as the data. The time resolution used in the model corresponds to the time resolution of the precipitation data and to the annual scale for the ice core data (i.e. yearly averages based on daily precipitation isotopic composition weighted by the amount of daily precipitation). The data type is identified as 1 for precipitation samples and 2 for ice core data.

Figure 8 Average seasonal cycles from precipitation data over the available period (dashed lines with points) and simulated by the ECHAM5-wiso model over the period 1979–2013 (solid lines) of the temperature (in ^∘C) (a), the precipitation (in mm w.e. yr⁻¹), the precipitation δ¹⁸O (in ‰) (c) and the deuterium excess (in ‰) (d). Data are shown for different durations depending on sampling, while model results are shown for the period 1979–2013. The number of points used for the observations is given in Table 3.

The simulated mean seasonal δ¹⁸O amplitude increases gradually from coastal regions to central Antarctica (more than 15 and up to 25 ‰ for some areas; Figs. 7a and 8c, solid lines). The model–data comparison suggests that this pattern is correct and that the model may underestimate the inland seasonal amplitude. As previously reported for annual mean values, systematic offsets are also identified for seasonal variations, with a systematic overestimation of monthly isotopic levels inland (e.g. for Dome C and Dome F) and a systematic underestimation on the coast (e.g. for Vernadsky and Halley). The model–data mismatch is largest during local winter months.

Minima are observed and simulated in winter (May–September) at most locations, except for Rothera and Vernadsky where the data show a minimum in July but the model produces a minimum in late autumn (April). Maximum values are observed and simulated in local summer (December–January); a secondary maximum is also sometimes observed and simulated in late winter (August–September). Data from Marsh station show maxima in January, April, and August, whereas the model only produces a single summer maximum value.

In summary, we report no systematic bias of the seasonal temperature amplitude (Fig. 8a). The seasonal pattern for the temperature is similar compared to δ¹⁸O, with minima in winter and the largest model–data mismatch in winter. Secondary minima or maxima cannot be discussed with confidence, as they have low amplitudes. We also highlight that model–data offsets are larger in winter. Note that precipitation and d seasonal cycles are described in Sect. 3.2.4.

3.2.3 δ¹⁸O–T relationships

Table 5 reports the temporal δ¹⁸O–T relationships established from precipitation and temperature observations and those simulated by ECHAM5-wiso. This calculation is based on daily or monthly values (depending on the sampling resolution) and includes seasonal variations. The data display significant linear relationships for all sites but Marsh (p value = 0.07), with an increased strength of the correlation coefficient from the coast (e.g. r = 0.38 at Rothera) to the East Antarctic Plateau (e.g. r = 0.88 at Dome F). The lowest slopes are identified in the peninsula region, with a mean slope of 0.32 ‰ ^∘C⁻¹ for Rothera and Vernadsky, while the highest slopes occur over the East Antarctic Plateau, with a mean slope of 0.68 ‰ ^∘C⁻¹ for Dome C and Dome F. These temporal slopes appear mostly lower than the spatial slopes and those expected from a Rayleigh distillation with a single moisture source (typically 0.8 ‰ ^∘C⁻¹).

Table 5 Slope (in ‰ ^∘C⁻¹), correlation coefficient, and p value of the δ¹⁸O–temperature linear relationship from precipitation measurements over the available period at daily or monthly (when the name of the station is associated with an asterisk) scale, depending on the time resolution of the data, and from the ECHAM5-wiso model over the observed period at the time resolution of the data. Numbers in brackets correspond to the standard errors.

In the ECHAM5-wiso model, as for the data, the simulated isotope–temperature relationship is statistically significant for all sites but Marsh (p value = 0.06). However, correlation coefficients are very small for Rothera and Vernadsky, which are thus excluded from further analyses. In the simulation, correlation coefficients are the highest for Halley, Dome C, and Dome F (up to 0.55) and the lowest for Neumayer (as low as 0.29). The slope is the lowest at Neumayer, with a value of 0.29 ‰ ^∘C⁻¹, increases at Halley with a value of 0.48‰ ^∘C⁻¹, and is the highest over the plateau with values of 0.70 ‰ ^∘C⁻¹ at Dome C and up to 0.94 ‰ ^∘C⁻¹ at Dome F.

To summarise, ECHAM5-wiso tends to underestimate the strength of the isotope–temperature relationship, but correctly simulates a larger strength of the correlation in the central Antarctic Plateau compared to coastal regions. There are significant differences in the isotope–temperature slopes for both coastal and central plateau locations. While there is some agreement (e.g. for Dome F and Halley), the model also produces non-realistic slopes, with a much larger slope than observed at Dome C, for instance.

3.2.4 The δD–δ¹⁸O relationship and d patterns

The δ¹⁸O–δD linear relationship is expected to be affected by different kinetic fractionation processes, for instance those associated with changes in evaporation conditions. We first compare the δ¹⁸O–δD linear relationship in the available precipitation and ice core data and simulated by ECHAM5-wiso (Table 6). Significant correlation is observed for all observational datasets but Marsh, as expected from meteoric samples, assuming correct preservation of samples and accurate isotopic measurements. We stress that the smallest correlation coefficient is identified at Vernadsky (r = 0.96), suggesting potential artefacts for this record. In the observations, the δD–δ¹⁸O slope varies across regions. While slopes higher than for the global meteoric waterline (i.e. > 8 ‰ ‰⁻¹) are identified at DDU and in Dronning Maud Land, lower slopes are identified in the Antarctic Peninsula (6.6 to 7.0 ‰ ‰ ⁻¹) and in the central East Antarctic Plateau (6.5 and 6.4 ‰ ‰⁻¹ at Dome C and Dome F, respectively). In the model, outputs also display significant linear relationships. They show higher values of the slope than observed in the Antarctic Peninsula, at DDU, and at Dome F and lower than observed for the other regions, including Dome C. These results appear coherent with associated coastal versus inland temperature and isotopic distillation biases.

Table 6 Slope (in ‰ ‰⁻¹), correlation coefficient, and p value of the δ¹⁸O–δD linear relationship from precipitation measurements (top of the table) and ice core data (bottom of the table) over the available period at daily or monthly scale (identified with an asterisk) and from the ECHAM5-wiso model over the observed period at the time resolution of the data for the precipitation and at the annual scale for the ice core data. Numbers in brackets correspond to the standard errors.

Figure 6 compares the spatial patterns of the d time-averaged model–data difference (characterised at 293 grid cells in our database; see Fig. 6a), and the situation is contrasted with 50 % of positive and negative differences. We can identify systematic trends, with an underestimation of the mean d levels in ECHAM5-wiso for the central East Antarctic Plateau and the peninsula and an overestimation above Victoria Land (Fig. 6a). Due to the temperature dependency of equilibrium fractionation coefficients leading to a gradual deviation from the meteoric waterline (calculated at the global scale, at which a coefficient of 8 results from the average equilibrium fractionation coefficients), d increases when temperature decreases (Masson-Delmotte et al., 2008; Touzeau et al., 2016). For central Antarctica, the d bias is thus consistent with the warm bias and the lack of isotopic depletion. The upper and lower quartiles of the model–data differences range within ±1.5 ± 0.1 ‰, suggesting that the model outputs remain close to those observed.

The d pattern is similar to that of δ¹⁸O: ECHAM5-wiso underestimates the d standard deviation for 90 % of grid cells, with an interquartile range comparable to the one for the ratio of standard deviations for δ¹⁸O (Fig. 6b). Table 7 displays the comparison of the statistics between d in the observations and in ECHAM5-wiso. In the observations, the time-averaged d is particularly low in the peninsula (−3.6 to 8.6 ‰), intermediate in the coastal regions of Dronning Maud Land, Victoria Land, and Adélie Land (4.4 to 8.6 ‰), and very high in the central Antarctic Plateau (up to 17.5 ‰ for Dome C). Lower coastal values and higher inland values are captured by ECHAM5-wiso, albeit with large offsets for each site reaching several per mille. These findings are consistent with the map showing the time-averaged precipitation d simulated by ECHAM5-wiso over the period 1979–2013 (Fig. 6a), with very low coastal values (close to zero) and increasing values towards the interior of Antarctica, reaching values higher than 16 ‰ on the plateau. ECHAM5-wiso mainly underestimates the d intra-annual standard deviation for 10 sites out of 15 (Table 7 and Fig. 6b).

Table 7 Mean value ± standard deviation (in ‰) of sub-annual d in observational time series at daily or monthly scale (identified with an asterisk) for the precipitation and for the ice core data and simulated d by ECHAM5-wiso for the same time period as the observations for precipitation and at the annual scale for the ice core. Mean values which are overestimated by ECHAM5-wiso are written in italic.

Figure 8d depicts the mean d seasonal patterns of the precipitation data and corresponding model outputs. The data show different patterns from one location to another. While d measured at Neumayer, Halley, and Rothera displays a maximum in autumn (March–April), it appears in late autumn (May) at Marsh and in winter (June–August) at Vernadsky. Maxima for central stations are observed later, in May–July for Dome C and July–September for Dome F. In short, most coastal areas are associated with a maximum d in autumn, while central areas are associated with a later maximum d, i.e. in winter or late winter, that is thus in anti-phase with δ¹⁸O and temperature. The seasonal amplitude increases from the coast to the plateau. In the model, for central areas, a first d maximum is simulated earlier than observed (February–March for Dome F and May–June for Dome C), followed by a second maximum in late winter (August for Dome F and September for Dome C). For coastal areas, the amplitude of the simulated d signal is too small to unequivocally estimate the timing of the maximum. Note the very low value simulated at DDU in July, which appears to be an outlier when comparing this value with the average modelled d value for all days in August 1973 (+5.9 ‰). No link emerges between the modelled seasonal patterns in d and in temperature (Fig. 8a), accumulation (Fig. 8b), or δ¹⁸O (Fig. 8c).

Table 8 The d mean seasonal amplitude (in ‰) calculated for precipitation at daily or monthly scale (identified with an asterisk) and sub-annual ice core data, as well as simulated by ECHAM5-wiso for the same time period as each record. The data type is identified as 1 for precipitation samples and 2 for ice core records. Amplitude values that are overestimated by ECHAM5-wiso are written in italic.

Finally, Table 8 reports the d mean seasonal amplitude values for the precipitation data and ice core records, as well as for the model outputs covering the observation. They clearly show an increase in d seasonal amplitude from the coast to the plateau (see also Fig. 7b), with values varying from 6.7 at Halley to 41 ‰ at Dome C. ECHAM5-wiso systematically underestimates the d mean seasonal amplitude when compared with precipitation data, while it systematically overestimates it when compared with ice core data (from 9.4 to 15.5 ‰), with the exception of the GIP ice core. Again, we cannot rule out a loss of amplitude in ice core data compared to the initial precipitation signal due to the temporal resolution and post-deposition effects.

3.3 Strength and limitations of the ECHAM5-wiso model outputs

The isotopic model–data time-averaged biases appear coherent with temperature. A warm bias over central East Antarctica and a cold bias over coastal regions lead to a too-low and too-strong isotopic depletion, respectively. Temperature and distillation biases also explain the underestimation of d above the central East Antarctic Plateau.

However, some characteristics are not explained by model skills for temperature. At sub-annual timescales, ECHAM5-wiso always overestimates the standard deviation of δ¹⁸O in precipitation (Table 3), but results for d are mixed (Table 7). ECHAM5-wiso always underestimates the seasonal amplitude of δ¹⁸O and d in precipitation but always overestimates the seasonal amplitude of δ¹⁸O and d in firn–ice cores (Tables 4 and 8). Differences between the model and firn–core data are at least partially due to diffusion processes, but no clear reason can be given for the other isotopic biases.

We do not find any clear link between other model biases for d and those for temperature or δ¹⁸O.

Sampling Antarctic snowfall remains challenging (Fujita and Abe, 2006; Landais et al., 2012; Schlosser et al., 2016; Stenni et al., 2016). Sampling is likely to fail to capture small events and may also collect surface snow transported by winds or hoar. Snow samples may undergo sublimation before collection. The fact that ECHAM5-wiso appears to overestimate the variability of precipitation isotopic composition may be related to an improper characterisation of the full day-to-day variability of real-world precipitation from daily precipitation sampling. Alternatively, this feature may also arise from a lack of representation of small-scale processes (boundary layer processes, wind characteristics, snow–atmosphere interplays) in ECHAM5-wiso. These processes may contribute to a local source of Antarctic moisture (through local recycling), reducing the influence of large-scale moisture transport (resolved by ECHAM5-wiso nudged to reanalyses) on the isotopic composition of precipitation and its day-to-day variability.

Caveats also limit the interpretation of the comparison of ECHAM5-wiso precipitation outputs with surface snow or shallow ice core data. Such records are potentially affected by post-deposition processes, such as wind scoring, erosion, snow metamorphism between precipitation events, and diffusion.

Our apparently contradictory findings for model–data comparisons with respect to inter-annual variations (from ice cores) and inter-daily variations (from precipitation data) call for more systematic comparisons between δ¹⁸O records of precipitation and ice cores at the same locations over several years.

4 Use of ECHAM5-wiso outputs for the interpretation of ice core records

Figure 9 Linear analysis of annual ECHAM5-wiso outputs from 1979–2013 for the temporal δ¹⁸O–temperature relationship (using the 2 m temperature and the precipitation-weighted δ¹⁸O). Maps show the slope of the linear regression (‰ ^∘C⁻¹) at the right side (a, c, e) and the correlation coefficient at the left side (b, d, f). The upper plots use outputs at the spatial scale (a, b), the middle plots at the inter-annual scale (c, d), and the lower plots at the seasonal scale (e, f). Areas where the results of the linear analysis are not significant are hatched (p value > 0.05).

In this section, we use the model outputs to help in the interpretation of ice core data: we quantify the inter-annual isotope–temperature relationships (Sect. 4.2) and characterise the spatial distribution of seasonal δ¹⁸O–d phase lag. Based on the confidence we can have in the model for each of the seven aforementioned regions (see Sect. 1 and Fig.1), we formulate recommendations for the future use of ECHAM5-wiso outputs (Sect. 4.3).

4.1 Spatial and temporal isotope–temperature relationships

First, we use ECHAM5-wiso to investigate spatial δ¹⁸O–temperature relationships (Fig. 9a and b) and then inter-annual (Fig. 9c and d) and seasonal relationships (Fig. 9e and f). For spatial relationships, the strength of the linear correlation coefficient is higher than 0.8. The spatial slope shows regional differences. It is generally smaller near the coasts (less than 0.8 ‰ ^∘C⁻¹), with the exception of Dronning Maud Land, and increases at elevations higher than 2500 m a.s.l., with values above 1.2 ‰ ^∘C⁻¹ in large areas. Furthermore, ECHAM5-wiso simulates the spatial heterogeneity of the gradient in the central East Antarctic Plateau around Dome C, Dome A, and Dome F. Such variability may arise from the simulated intermittency of precipitation and from differences in condensation versus surface temperature.

At the inter-annual scale (Fig. 9c and d), results are not significant for large areas encompassing the Dronning Maud Land region, the Antarctic Peninsula, the Transantarctic Mountain region, the Ronne and Filchner ice shelf regions, part of Victoria Land, and along the Wilkes Land coast. For the whole continent, the correlation coefficient varies between 0.5 and 0.6 (with few values reaching 0.6 at the upper limit and 0.3 at the lower limit). Where correlations are significant, the inter-annual δ¹⁸O–temperature slope increases from the coasts (0.3 to 0.6 ‰ ^∘C⁻¹) to the inland regions, where it can exceed 1 ‰ ^∘C⁻¹ for some high-elevation locations. The low correlation may be due to the small range of mean annual temperature over the period 1979–2013 and is not necessarily indicative of a weak sensitivity to temperature change.

Finally, at the seasonal scale, results are significant almost over the whole continent (with the exception of two little areas in the peninsula and East Antarctica) and the correlation coefficients are equal to 1 everywhere but along the coastal regions in the Indian Ocean sector, where the correlation coefficient can decrease down to 0.75. Slopes are lower than for spatial and inter-annual relationships, with values from 0.0 to 0.3 ‰ ^∘C⁻¹ along the coast (higher over Dronning Maud Land and the Ross Ice Shelf region), around 0.5 ‰ ^∘C⁻¹ inland for altitudes lower than 2500 m a.s.l. (with the exception of lower values above the Transantarctic Mountains), and up to 0.8 ‰ ^∘C⁻¹ over the East Antarctic Plateau.

To conclude, the coherent framework provided by the ECHAM5-wiso simulation covering the period 1979–2013 shows that annual δ¹⁸O and surface temperatures are only weakly linearly related in several areas. This suggests that the inter-annual variability of δ¹⁸O is controlled by other processes, for instance those associated with synoptic variability and changes in moisture source characteristics (Sturm et al., 2010; Steiger et al., 2017). Moreover, our results rule out the application of a single isotope–temperature slope for all Antarctic ice core records on the inter-annual timescale, and the seasonal isotope–temperature slope is not a surrogate for scaling inter-annual δ¹⁸O to temperature.

We have also used the simulation to explore linear relationships between d and surface air temperature, without any significant results (not shown).

4.2 δ¹⁸O–d phase lag

Deuterium excess (d) has originally been interpreted as a proxy for relative humidity at the moisture source (Jouzel et al., 2013; Pfahl and Sodemann, 2014; Kurita et al., 2016). However, recent studies of Antarctic precipitation data combined with back-trajectory analyses did not support this interpretation (e.g. Dittmann et al., 2016; Schlosser et al., 2017), calling for further work to understand the drivers of seasonal d variations. The phase lag between d and δ¹⁸O was initially explored to identify changes in evaporation conditions (Ciais et al., 1995). In ECHAM5-wiso, this phase lag is calculated as the lag that gives the highest correlation coefficient between d and δ¹⁸O (Fig. 10) using the mean seasonal cycle from monthly averaged values. If there were no seasonal change in moisture origin and climatic conditions during the initial evaporation process, one would expect d to be in anti-phase with δ¹⁸O due to the impact of condensation temperature on equilibrium fractionation. For regions with small seasonal amplitude in condensation temperature, a constant initial isotopic composition at the moisture source would imply a stable d year-round. In such regions, the simulated phase lag likely therefore reflects seasonal changes in the d of the initial moisture source. The comparison with precipitation data (Sect. 3.2.4) showed that ECHAM5-wiso had low seasonal amplitude in coastal regions (Fig. 8d), making the discussion of seasonal maxima difficult. These comparisons are also limited by the duration of the precipitation records. Here, we use the full simulation (1979–2013) to investigate the phase lag between the mean seasonal cycle of d and δ¹⁸O. Clear spatial patterns are identified for the distribution of this phase lag (see Fig. 10). At intermediate elevations (between 1000 and 3000 m a.s.l.), d seasonal variations occur in phase (within 2 months) with the seasonal cycle of δ¹⁸O (and surface air temperature). By contrast, a phase lag of several months is identified over coastal areas and over the central East Antarctic Plateau. Along the Wilkes Land coast and the Dronning Maud Land region, the time lag is between 2 and 4 months below 1000 and 500 m a.s.l., respectively. Over the West Antarctic Ice Sheet, the phase lag is higher than 2 months below 500 m a.s.l. and can even reach 6 months (indicating an anti-phase between d and δ¹⁸O). Over the central East Antarctic Plateau (above 3000 m of elevation), the phase lag reaches several months again, especially near Dome C. Obtaining longer precipitation records and comparing the phase lag identified in precipitation and surface snow records would be helpful to understand whether post-deposition processes, which are not included in ECHAM5-wiso, affect this phase lag. The different characteristics of seasonal d changes suggest different seasonal changes in moisture origin at coastal, intermediate, and central plateau regions, supporting the identification of specific coastal versus inland regions to assess the isotope–temperature relationships. Note that the few available datasets are in line with the simulation.

Figure 10 Best correlated phase lag between the mean seasonal cycle of deuterium excess and that of δ¹⁸O simulated by ECHAM5-wiso over the period 1979–2013 (colour shading) and calculated from precipitation data (circles). The sign provides information on the sign of the correlation between δ¹⁸O and d (e.g. positive numbers correspond to a correlation, while negative numbers correspond to an anti-correlation). The absolute value corresponds to the lag (in months) between δ¹⁸O and deuterium excess corresponding to the highest correlation of monthly averaged values. This figure also displays the Antarctic topography, with isohypses (in m a.s.l.).

4.3 Recommendations for the different regions of Antarctica

In this section, we summarise our findings based on the model–data comparisons and the analysis of model outputs for the seven Antarctic regions selected by the Antarctica2k program, as shown in Fig. 1 (Stenni et al., 2017). The regions depend on geographical and climatic characteristics. Results from Sect. 3 were averaged over each region and are given in Table 9.

Table 9 Evaluation of the ECHAM5-wiso model for seven Antarctic regions: East Antarctic Plateau, coastal Indian, Weddell Sea, peninsula, West Antarctic Ice Sheet, Victoria Land, and Dronning Maud Land (7). We regionally averaged the time-averaged δ¹⁸O mean (model–data) differences (in ‰), the inter-annual δ¹⁸O standard deviation (model ∕ data) ratio, the time-averaged d mean (model–data) differences (in ‰), and the inter-annual d (model–data) standard deviation ratio using only precipitation data. Italic cells correspond to parameters for which we support the validity of the use of ECHAM5-wiso for the considered region, underlined cells to parameters for which we suggest some caution, and bold cells to parameters for which we suggest not using ECHAM5-wiso outputs for the considered region. The numbers in brackets correspond to the number of data points.

We first discuss the systematic model biases. The maximum time-averaged model–data differences (3.8 and 2.6 ‰ for δ¹⁸O and d, respectively) are identified in the Weddell Sea area. Minimum time-averaged model–data differences occur in different regions for δ¹⁸O and d (Victoria Land and Dronning Maud Land, respectively).

For inter-annual standard deviation, the model–data mismatch is smallest for Victoria Land (ratio of 1.1 and 1.0 for δ¹⁸O and d, respectively). Results for δ¹⁸O show that the simulated inter-annual variability can be considered close to reality (model–data ratio higher than 0.7) only for Victoria Land and the plateau, acceptable (model–data ratio higher than 0.5) for the Weddell Sea area and the West Antarctic Ice Sheet, but significantly different from observations in the other three regions. The model–data mismatch is larger for d inter-annual variability, with acceptable inter-annual variability only for Victoria Land and the plateau. However, these results are clearly limited by the low number of observational records for some regions.

Table 10 provides a brief overview of ECHAM5-wiso outputs for our seven regions of interest in terms of mean climate and isotopic variables, their standard deviation, seasonal amplitude, and the calculated regional δ¹⁸O–T relationship. The main findings are again the highest slope simulated for the central Antarctic Plateau, followed by the Dronning Maud Land and West Antarctic Ice Sheet regions, and weak correlations in some regions (Weddell Sea, Antarctic Peninsula) where water stable isotope outputs are not good predictors of inter-annual temperature change within ECHAM5-wiso, together with low correlations and slopes for the other coastal regions (Indian Ocean sector, Victoria Land).

Table 10 Exploration of the ECHAM5-wiso model outputs (1979–2013) for seven Antarctic regions: east plateau, coastal Indian, Weddell Sea, peninsula, West Antarctic Ice Sheet, Victoria Land, and Dronning Maud Land (7). For each of the variables precipitation (in mm w.e. yr⁻¹), temperature (in ^∘C), δ¹⁸O (in ‰), and d (in ‰), we regionally averaged the annual mean values (lines 1 to 4), the inter-annual standard deviation (lines 5 to 8), and the mean seasonal amplitude (lines 9 to 12). Finally, we calculated the statistics of the inter-annual δ¹⁸O–temperature linear relationship: the slope (noted as “a”, in ‰ ^∘C⁻¹), the correlation coefficient (noted as “r”), and the p value for each region.

5 Conclusions and perspectives

This study presents a systematic evaluation of a present-day Antarctic climate simulation using the ECHAM5-wiso atmospheric circulation model equipped with water stable isotopes. For this simulation covering the period 1960–2013, the model has been nudged to ERA atmospheric reanalyses. In particular, we tested its ability to correctly capture time-averaged values, inter-annual variations, and seasonal cycles in surface mass balance, temperature, and precipitation isotopic composition in Antarctica. As fare as possible, we discarded model results prior to 1979, as model–data differences prior to 1979 may arise from uncertainties in the reanalyses prior to the period for which satellite data were assimilated.

Despite some divergences, simulated P-E values are found to be a good surrogate for SMB. Most artefacts in modelled δ¹⁸O are coherent with those for temperature, with systematic biases in different regions. Some of these artefacts may be linked to the nudging method and the reanalyses. Model–data comparisons are limited by data availability and by the fact that deposition and post-deposition processes are not considered in the simulation. This is particularly true for precipitation amounts, for which there is a lack of direct measurements, and isotopic analysis for many regions at a multi-annual timescale. A systematic comparison between water isotope measurements from precipitation and surface snow or ice core samples is needed for further in-depth studies of this topic. We note a lower quantitative performance from ECHAM5-wiso for d (time-averaged values and inter-annual standard deviations) than for δ¹⁸O, beyond its remarkable ability to resolve the spatial distribution of time-averaged d values. Our findings confirm several other studies conducted in other regions highlighting the fact that atmospheric models including ECHAM5-wiso tend to underestimate the variability of d in surface vapour (e.g. Steen-Larsen et al., 2016). Expanding earlier site-specific studies, we show that the strength and slope of the δ¹⁸O–temperature linear relationship is dependent on the timescale in Antarctica over the last 4 decades. This finding has implications for past temperature reconstructions using ice core records. Finally, interesting results emerge for regional differences in the phase lag between the mean seasonal cycle in δ¹⁸O and d, calling for further studies to better characterise this feature in precipitation and ice core records and better understand the implications of these lags for the representation of seasonal changes in moisture source effects.

Our study deserves to be expanded to other atmospheric models equipped with water stable isotopes and other nudged simulations using different reanalyses datasets to assess the robustness of our findings. Furthermore, obtaining more high-resolution ice core records is crucial to better assessing model skills for inter-annual variations. More measurements of precipitation, surface snow, and vapour monitoring for water isotopes would also help to better characterise deposition and post-deposition processes, their implication for model–data evaluation studies, and for an improved climatic interpretation of ice core records.