Air mass origin and local impacts on Antarctic snow isotopic composition: an observation and modelling study

Water stable isotopes (δ¹⁸O and δD) from ice cores are widely used to reconstruct past temperature variations through their well-established relationship with local air temperature, commonly referred to as “isotopic paleothermometer”. However, depositional and post-depositional effects lead to large uncertainties in the use this proxy in Antarctica. The magnitude of these uncertainties strongly depends on site location, with larger impacts in low-accumulation regions of East Antarctic Plateau. Depositional effects include origin of moisture, which exhibits asymmetries shaped by the continent's geographical and topographical features, as well as precipitation intermittency, which introduces aliasing in the archived signal. Post-depositional processes, such as sublimation and firn-atmosphere exchange, further alter the isotopic composition of snow before its transformation into ice, thereby modifying the correlation between δ¹⁸O and temperature. Here, we present new water isotope measurements from surface snow collected during the East Antarctic International Ice Sheet Traverse (EAIIST), across a remote region of the East Antarctic Plateau. The traverse – crossing a transitional zone between predominately Indian and Pacific moisture sources – provides unique insights into the key role of air mass origin in shaping the isotopic composition of snow. Comparison with LMDZ6iso simulations indicates that the model successfully captures the spatial variability of δ¹⁸O-temperature relationship between different basins, with statistically significant correlations (p < 0.05) when the analysis is extended to the Antarctic dataset. This agreement further suggests the model's ability to predict the temporal slope required for calibrating isotopic ice-core records used for temperature reconstructions, even in regions influenced by multiple moisture sources. Temporal slopes based on monthly precipitation values range from 0.4 ‰ °C⁻¹ to 0.5 ‰ °C⁻¹ for the EAIIST drilling sites. Finally, we quantify the impact of sublimation on isotopic composition of surface snow. Including sublimation in the modelling of surface snow reduces the discrepancy between observed and modelled values, compared to simulations accounting precipitation, from 1.9 ‰ to 1.3 ‰ for δ¹⁸O and from 6.6 ‰ to 2.9 ‰ for d-excess. These results highlighting the key role of this post-depositional process on the Antarctic Plateau.

1 Introduction

Ice cores are valuable archives of past climatic and environmental conditions through many environmental proxies trapped in the ice, such as air bubbles and aerosols, as well as the ice itself, such as its stable water isotopic composition. In East Antarctica, water stable isotopes (δ¹⁸O and δD) from the deepest ice cores have been used to reconstruct past temperatures extending back 800 000 years (EPICA community members et al., 2004; Jouzel et al., 2007), and recent efforts aim to extend these reconstructions to 1.5 million years (Beyond EPICA Oldest Ice Core – Parrenin et al., 2017; (Lilien et al., 2021). In this region, ice core temporal resolution is typically limited to multi-annual to multi-decadal scale (Ekaykin et al., 2002; Baroni et al., 2011; Münch et al., 2016). This limitation is due to low accumulation rates and intermittency of precipitation which cause aliasing of the climatic signal recorded by sporadic snowfall events (Ekaykin et al., 2016; Laepple et al., 2018; Casado et al., 2020; Münch et al., 2021). However, ice cores in high-accumulation coastal regions of Antarctica cover shorter temporal ranges but provide finer temporal resolution that can capture seasonal-scale variations (Casado et al., 2023).

Several studies have reported the empirical linear relationship between the isotopic composition of snow and local temperature (Dansgaard, 1964; Lorius et al., 1969; Touzeau et al., 2016). Nevertheless, this relationship is not stable in time and space, as widely documented by direct snow observations and model results (Goursaud et al., 2018; Jouzel et al., 1997). These variations are related to changes in evaporative conditions and transport pathways to condensation sites which influence the isotopic equilibrium and kinetic fractionation processes at each step of the distillation trajectories (Charles et al., 1994; Masson-Delmotte et al., 2011; Werner et al., 2011, Casado et al., 2017). In addition, the second order parameter deuterium excess (d-excess = δD − 8 × δ¹⁸O, (Dansgaard, 1964) is highly sensitive to kinetic effects occurring both during evaporation at the ocean surface and along distillation pathways during atmospheric transport. As a result, while over the ocean d-excess is very sensitive to near-surface relative humidity and weakly to temperature, it exhibits the opposite behaviour in central East Antarctica, where it may reflect, at least partially, different moisture origins (Petit et al., 1991; Risi et al., 2013; Vimeux et al., 1999; Jouzel et al., 2003; Masson-Delmotte et al., 2005; Neumann et al., 2005).

Additionally, the reliability of water isotopes as climate proxy from ice core is debated due to the local influence of post-depositional processes in the snowpack. These processes can alter the snow isotopic composition during prolonged exposure at the atmosphere-snow interface and modify the original precipitation's isotope-temperature relationship (Petit et al., 1982; Stenni et al., 2016, Touzeau et al., 2016, Casado et al., 2018). Key post-depositional mechanisms include wind-driven snow redistribution and vapor exchange with the atmosphere, such as sublimation-condensation processes, snowmelt and vapor diffusion within the snowpack driven by forced ventilation (Steen-Larsen et al., 2014; Beria et al., 2018; Casado et al., 2018; Wahl et al., 2022; Ollivier et al., 2025). Among these, sublimation plays an important role especially during the warmer months, decreasing d-excess values (Landais et al, 2017, Casado et al, 2021).

These post-depositional alterations are poorly understood and often overlooked in climatic reconstructions, underscoring the need to better constrain the processes shaping the isotope-temperature relationship in surface snow to improve the reliability of ice core-based temperature estimates (Xiao et al., 2013; Ma et al., 2020).

In this study, we investigate the isotopic composition of snow samples collected roughly every 20 km during the East Antarctic International Ice Sheet Traverse (EAIIST – (Traversa et al., 2023) in 2019–2020. The dataset includes both surface snow and bulk samples, which reflect precipitation over time scales ranging from seasonal to multi-annual periods, respectively. This scientific traverse extended from the coast in Adélie Land into the high interior of the East Antarctic Plateau. We assess the influence of Indian vs. Pacific Ocean moisture sources on the spatial variability of the relationship between δ¹⁸O composition of snow and the 2 m air temperature. In addition, we present the d-excess vs. δ¹⁸O relationship, which reveals different supersaturation pathways under cold conditions. To place our results in a broader context, we compare them with a subset of surface snow samples from Indian and Pacific sectors of Antarctica, drawn from the dataset compiled by (Masson-Delmotte et al. (2008).

We then evaluate the ability of the atmospheric general circulation model LMDZ6iso to reproduce the observed δ¹⁸O-temperature spatial slope by comparing surface snow observations with modelled precipitation. Subsequent to the spatial analysis, we further assess the model's capability to reproduce temporal slopes, which are commonly used in paleoclimate reconstructions of temperature from ice cores.

Finally, surface snow samples collected along the plateau during the outward and return paths of the traverse provide an opportunity to assess post-depositional isotopic modification over time. By comparing the observations with precipitation model outputs, we isolate the impact of sublimation effect on the snow during summer months. This quantification is achieved by applying a metamorphism model proposed by Casado et al. (2021).

Figure 1 (a) Map of EAIIST traverse and (b) precipitation period derived from ERA5, corresponding to surface and bulk snow samples across the traverse. Main sampling sites are marked with red squares: Dumont D'Urville (DDU), Dome C (DC), AGO5, PALEO, and Megadune (MD).

2 Methods

2.1 Geographical area EAIIST

The EAIIST traverse took place during the austral summer 2019–2020. The traverse started in coastal Adélie Land, near the high accumulation site of Dumont D'Urville (DDU: 66°39^′ S, 140°01^′ E; elevation = 47 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$; mean annual temperature of −10.8 °C), until the interior plateau site of Megadune (MD: 80°35^′ S, 121°35^′ E; elevation = 2973 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$), and then returned along the same route back to DDU (Fig. 1a). The traverse route passed by the Dome C site (DC: 75°06^′ S, 123°23^′ E; elevation = 3233 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$; mean annual temperature of −54.5 °C) where the EPICA ice core was drilled (EPICA community members et al., 2004) and where the Concordia station is currently located. Dome C represents the maximum elevation of the traverse.

2.2 EAIIST snow samples and isotope measurements

Two types of surface snow samples were collected: 85 surface samples, representing the upper 3 cm of snow, and 52 bulk samples, consisting of snow integrated over a vertically dug 1 m-deep snowpit. Surface samples were taken at each stop during daytime approximately every 20 km along the 1600 km route from DDU to MD, and onward to DC. A single surface sample at one location was obtained by mixing snow from 10 locations randomly selected over a similar area of 100 m × 100 m, representing an average. Bulk sampling required longer processing time digging 1 m deep snowpit and scratching the snow along the vertical profile and was therefore carried out only during lunch and evening stops. The samples were collected with 50 mL Corning tubes which were sealed to prevent air exchange and kept frozen until arrival to laboratories in Europe. Samples were distributed between the Laboratoire des Sciences du Climat et de l'Environnement of Paris (LSCE, France) and Ca' Foscari University of Venice (UNIVE, Italy) for water isotopic analysis. The samples were melted only immediately prior the measurement to minimise potential alterations of the water isotopic composition. Analyses were performed using a Cavity-Ring Down Spectroscopy (CRDS) analyser PICARRO model L2130-i at UNIVE and LSCE.

The isotopic composition of snow is expressed in delta-notation (‰) (Craig, 1961) relative to laboratory standards, which were previous calibrated against the international standards V-SMOW (Vienna Standard Mean Ocean Water) and V-SLAP (Vienna Standard Light Antarctic Precipitation). The accuracy of PICARRO measurements was determined as the mean difference between measured and true values of laboratory standards, with uncertainty represented by their standard deviation. This yielding an accuracy of −0.01 ‰ for δ¹⁸O, −0.07 ‰ for δD, and −0.02 ‰ for d-excess, with corresponding uncertainties of ± 0.07 ‰, ± 0.4 ‰, and ± 0.4 ‰. Previous inter-calibration experiments revealed mean discrepancy between UNIVE and LSCE measurements of the same samples equal to 0.14 ‰ and 0.80 ‰, for δ¹⁸O and d-excess, respectively (Petteni et al., 2025).

Figure 2 Precipitation-weighted air mass back-trajectory distributions for the period 2009–2019 at four key sites: DC, AGO5, PALEO, and MD. Colors indicate air mass origins: Indian (red), Pacific (blue) and Atlantic (orange) oceans.

To estimate the order of magnitude of spatial isotopic variability in surface snow at one location, we report the standard deviation (SD) from previous Antarctic datasets. At Dome C (Casado et al., 2018), the local variability is around 3.4 ‰ for δ¹⁸O and 4.1 ‰ for d-excess (defined as 2 SD of replicates obtained in an area of 100 m × 100 m). Since the surface samples correspond at mixed snow from 10 locations, the impact of spatial variability on the error of the mean is expressed as the standard deviation of the mean (also known as standard error, SE):

$$\begin{array}{}\text{(1)}& \text{SE}={\displaystyle \frac{s}{\sqrt{n}}}\end{array}$$

where s is the SD calculated from individual samples and n is the number of subsamples. Applying this, the expected spatial variability for our composite surface samples is $\mathrm{3.4}\mathrm{‰}/\sqrt{\left(\mathrm{10}\right)}$ ∼ 1.1 ‰ for δ¹⁸O and $\mathrm{4.1}\mathrm{‰}/\sqrt{\left(\mathrm{10}\right)}$ ∼ 1.3 ‰ for d-excess. Similar SE are observed in the coastal area of Dumont D'Urville, based on the top 2.5 cm of snowpit samples reported by (Landais et al., 2017). Overall, while these values can vary across sites separated by hundreds of kilometres, we set an indicative uncertainty arising from local spatial variability of 1.1 ‰ for δ¹⁸O and 1.3 ‰ for d-excess, as reference for interpreting our data.

2.3 Antarctic dataset

To provide a broader spatial analysis, we compare our data with the Antarctic surface snow database of Masson-Delmotte et al. (2008), which offers a comprehensive overview of isotopic variability across the continent. From the original dataset, which includes different types of snow samples, we selected surface snow, bulk snow, snowpit and firn cores that capture signal ranging from annual to approximately 20 years, based on sample depth and local precipitation rates. For the comparison with our snow samples, we divide the dataset in Pacific and Indian sectors, based on (Sodemann and Stohl, (2009). To do this, we classified the region at west of 60° W and the area near the Ross Sea as Pacific sector, while the sampling sites located north then 80° S and between 60 and 180° E as Indian sector. For this dataset, mean δ¹⁸O values from LMDZiso model were simulated over the 1980–2000 period, due to the lack of precise temporal coverage of the snow samples.

2.4 FLEXPART back-trajectories

The Lagragian particle dispersion model FLEXPART (FLEXible PARTicle) is employed to calculate 10 d back-trajectories of air masses at 12 h temporal resolution for each sampling site of the traverse located on the Antarctic Plateau. These trajectories are derived for the 500 hPa pressure level, representing mid-tropospheric transport pathways reaching the plateau.

With two trajectories calculated per day over the 2009–2019 period, we computed approximately 7300 trajectories for each site. To allow for a comprehensive multi-annual evaluation and emphasise the air masses leading precipitation events at the sampling sites, the trajectories were averaged and weighted by ERA5 precipitation rates (see Sect. 2.5). The classification of air mass origins into Indian, Pacific and Atlantic sectors is based on their latitudinal and longitudinal distributions, as illustrated in Sect. 3.1, Fig. 2.

2.5 ERA5 climatic signal

We use ERA5 reanalysis data (provided by the European Centre for Medium-Range Weather Forecasts) to estimate the precipitation interval and associated temperature conditions related to snow precipitation of surface and bulk snow samples (Fig. 1b). ERA5 provides hourly precipitation and temperature values at 0.25° spatial resolution.

Prior to the comparison with ERA5, all snow samples were converted to water equivalent using density values of trench measurements at Dome C (Ooms et al., 2026). The densities are equal to 290 kg m⁻³ for the upper 3 cm of snow and 320 kg m⁻³ for the upper 1 m. For each sample, precipitation events in ERA5 were sequentially accumulated until the target water-equivalent thickness of the sample was reached. Due to the strong gradient in precipitation rates from the coast (100–300 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$) to the plateau (20–50 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$), surface samples near DDU represent from 1 to 3 months of snowfalls, whereas those collected on the plateau correspond to up to ∼ 6 months. Similarly, bulk samples represent approximately 1 year of precipitation in coastal areas and up to 15 years at the highest-elevation sites. This estimate considers only snowfall-driven precipitation – excluding accumulation coming from clear-sky precipitation, such as diamond dust and vapor condensation (Stenni et al., 2016) or effects from wind redistribution, erosion, and sublimation – but provides a distinction between seasonal (surface) and multi-seasonal signals (bulk) signals recorded in the samples.

To investigate the δ¹⁸O-temperature (δ¹⁸O-T) relationship, we compare two temperature metrics. The first is the average 2 m temperature (T) which mainly captures the spatial climatic variability across the continent. The second is the precipitation-weighted temperature (T_pw) which has a more temporal emphasis reflecting the thermal conditions during snowfall events. Indeed, T_pw is an average temperature where the values are weighted according to the amount of snowfalls during each time interval. We calculate T and T_pw for each sampling site, based on the precipitation interval corresponding to surface snow samples (see Sect. 3.1). Furthermore, we calculate the average multi-annual (1980–2020) temperature for both EAIIST and Antarctic snow isotopic database (Masson-Delmotte et al., 2008), enabling a consistent spatial analysis for samples collected during different campaigns.

2.6 LMDZ6iso model data

LMDZ6iso (Risi et al., 2010) is the isotope-enabled version of the atmospheric general circulation model (AGCM) LMDZ6 (Hourdin et al., 2020). We use the LMDZ6iso version 20231022.trunk with the NPv6.3 physical package (Hourdin et al., 2023), which is nearly identical to the IPSL-CM6A atmospheric setup (Boucher et al., 2020), used for phase 6 of the Coupled Model Intercomparison Project CMIP6 (Eyring et al., 2016). We use LMDZ6's standard horizontal Low Resolution (LR) longitude-latitude grid (144 × 142), which corresponds to a 2.0° resolution in longitude and 1.67° in latitude. The vertical grid comprises 79 levels, with the lowest atmospheric level approximately 7 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ (above ground level) at Dome C. The simulation is nudged towards 6-hourly three-dimensional fields of temperature and wind from the ERA5 reanalysis (Hersbach et al., 2020), using a relaxation time scale of 3 h. Nudging is excluded below the sigma-pressure level equivalent to 850 hPa above sea level, allowing the physics and dynamics of the model to operate freely within the boundary layer. Surface ocean boundary conditions are derived from ERA5 reanalysis monthly mean sea surface temperature and sea-ice concentration fields. The simulation used in this study is described and evaluated over Antarctica in (Dutrievoz et al., 2025). Snow samples are simulated by stacking precipitation events until reaching the target thickness, corresponding to surface and bulk sample depths. The snow layers are precipitation-weighted before averaging, ensuring consistency when compared with the surface and bulk snow data. The error in LMDZ6iso precipitation at Concordia, defined as the difference between modelled and observed values for a supersaturation parameter λ = 0.004 K⁻¹, is 4.0 ‰ for δ¹⁸O and 1.9 ‰ for d-excess.

2.7 Snow metamorphism model

The snow metamorphism model proposed by Casado et al. (2021) describes the relative isotopic variations in snow induced by different processes. Additionally, the model quantifies the flux of water transferred from the snow to the atmosphere, starting from the isotopic composition of surface snow remaining after sublimation R_sub-snow, as described the following Eq. (2):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{R}_{\text{sub-snow}}^i=& {\mathit{\alpha }}_{\text{eq-sub}}^i\left(T\right)\\ & \times \left[R_{\text{snow}}^i{\left({\displaystyle \frac{D}{D^i}}\right)}^k\right(\mathrm{1}-\text{RH})+\text{RH}\cdot R_{\mathrm{a}}^i]\end{array}\end{array}$$

Where R_snow is the initial composition of the snow, Dⁱ and D are the kinetic diffusivity coeffients, α_eq-sub is the fractionation coefficient related to sublimation, RH is the relative humidity (Merlivat and Jouzel, 1979) and R_a is the isotopic composition of atmospheric vapor. The exponent k is the roughness parameter (Craig and Gordon, 1965), here set to 0.4. In particular, the model predicts that sublimation during the warmest months lead to a $\text{d-excess}/{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ slope steeper than −2 ‰ ‰⁻¹, whereas precipitation input follows a slope of −0.5 ‰ ‰⁻¹.

The isotopic difference between R_snow and R_{sub  snow} (expressed as ratio) is used to compute the flux of water transferred from the snow to the atmosphere. This computation takes into account the fractionation coefficient associated to sublimation and the initial water quantity of the snow sample. This approach will be used in Sect. 3.6, to isolate the impact of sublimation in surface snow, comparing the outbound and return samples of the traverse.

In addition, we simulate the final composition of the snow, starting from the outbound values and incorporating the snow precipitation input from LMDZ6iso model and taking into account the contribution of sublimation effects for the number of precipitation-free days.

3 Results

3.1 Air mass back-trajectories

We present precipitation-weighted air mass back-trajectories over the decadal period 2009–2019 for four representative plateau sites (Fig. 2). At Dome C, 90 % of the back-trajectories originates from the Indian Ocean (red), as expected for this part of the East Antarctic Plateau (Sodemann and Stohl, 2009). Although the traverse is topographically divided at Dome C, the transition between the predominance of Indian and Pacific air mass origins occurs farther south. At AGO5, approximately 70 % of precipitation events are associated with the Indian sector, with a secondary influence from Pacific Ocean (blue). Southern sites like PALEO and MD, are influenced by the Pacific Ocean, which accounts for over 85 % of precipitation events, with minor contributions from the Indian (5 %) and Atlantic (10 %) sectors. Based on these observations and the findings of Sodemann and Stohl, (2009) we set the transition between predominately Indian-sourced to predominately Pacific-sourced air masses at approximately 78° S along the EAIIST route (see Appendix A, Fig. A1). Note that these air-mass origins are more representative of the bulk samples rather than the surface snow, as they are computed as the average over 10 years.

3.2 ERA5 temperature signal

The δ¹⁸O of surface samples (top 3 cm) collected during the outbound and return ways of the traverse, which represent 5–6 months of snowfalls, exhibits strong correlation with both T and T_pw (R² = 0.94). However, the two relationships yield different slopes of 0.9 ‰ °C⁻¹ and 0.6 ‰ °C⁻¹, respectively (Fig. 3).

Figure 3 δ¹⁸O-T vs. δ¹⁸O-T_pw relationships evaluated for surface snow samples, calculating mean temperature and precipitation-weighted temperature over time intervals corresponding to 3 cm of snow precipitation.

The δ¹⁸O-T slope is consistent with the value obtained from surface, bulk and shallow snow samples of the Antarctic database (Masson-Delmotte et al., 2008) for Adélie Land based on multi-annual temperature (1980–2020). This similarity validates the use of multi-annual temperature (referred as T in the next sections) for spatial analysis of snow samples that represent at least an annual signal. On the other hand, the δ¹⁸O-T_pw slope is very close to that observed for precipitation at Dome C (0.5 ‰ °C⁻¹; Dreossi et al., 2024), suggesting that the δ¹⁸O-T_pw relationship is closer to approximating the temporal relationship between snow isotopes and local weather.

3.3 Spatial distribution of bulk isotopic composition

The snow isotopic composition of bulk samples ranges from −20.0 ‰ at the coastal site near DDU to −52.5 ‰ at the most inland Megadune site (a comparison with surface samples is provided in Appendix B, Fig. B1). This range greatly exceeds the uncertainty expected from spatial variability per site (∼ 1.1 ‰, see Sect. 2.2), and these values fall well within the range of previous Antarctic surface snow isotopic composition (Fig. 4). Between DDU and Dome C, temperature and δ¹⁸O both decline in parallel as the altitude increases. Yet, south of Dome C toward Megadune, the temperature increases as the altitude decreases, but δ¹⁸O slightly decreases. Similarly, the d-excess increases linearly from DDU to Dome C but then remains constant further south. Following our back-trajectory air mass results, we divide the EAIIST snow isotopic dataset into two sections split at 78° S.

Figure 4 Map of Antarctica showing the distribution of the δ¹⁸O values of EAIIST bulk samples for Section I (from DDU to DC) and Section II (from DC to MD), outlined with black contours. Values are shown alongside surface, bulk and shallow snow samples from the Antarctic dataset (Masson-Delmotte et al., 2008).

Figure 5 Spatial distribution of δ¹⁸O and d-excess values of bulk samples for Section I (red dots – north of 78° S) and Section II (blue dots – south of 78° S). The vertical bars represent the uncertainty associated with local variability. Temperature, elevation and distance from DDU are reported.

Table 1 Slope and correlation coefficients (R²) of linear relationships between isotopic composition and geographical/climatic variables, calculated for the Antarctic dataset and Sections I and II of the EAIIST traverse. Relationships are shown in bold when statistically significant (defined by p < 0.05).

Section I (north of 78° S) exhibits the expected linear relationship between isotopic composition of snow and local climatic and geographic parameters (Table 1, Fig. 5). These factors include distance from the nearest coast, latitude, elevation and mean annual 2 m air temperature. In contrast, in Section II (south of 78° S), these linear correlations disappear for all the mentioned parameters.

Within Section II, where no abrupt temperature variations are observed, two anomalous sites are notable: one located 20 km south of PALEO site, with a δ¹⁸O value of −47.9 ‰, and another at the Windcrust site (WC) with a value of −48.9 ‰, roughly 4 ‰ higher than those of the nearby samples. Both sites are also characterised by extremely low d-excess values of 7.0 ‰ and 3.3 ‰, respectively, compared to those observed on the plateau. Such values likely indicate local sublimation (Casado et al., 2021; Wahl et al., 2022). Supporting this interpretation, field observations indicate the presence of large hoar crystals at these sites which experience relatively large strong wind, potentially of katabatic origin (Bintanja, 1998). Although weaker on the plateau than along the continental margins, these winds can persists for extended periods, creating dry conditions that promote sublimation during summer (Grazioli et al., 2017). In Sect. 3.6, we further investigate the impact of sublimation on surface snow based on the comparison between the samples collected during the outward and return paths of the traverse.

Figure 6 (a) δ¹⁸O-T relationship and (b) d-excess vs. δ¹⁸O for bulk samples of the traverse compared with Antarctic database for Pacific and Indian sectors.

3.4 Moisture origin impact on isotopic composition: observations vs. model outputs

We present the EAIIST and Antarctic datasets subdivided into Pacific and Indian sectors. For the Antarctic dataset, the linear regressions for the δ¹⁸O-T relationship independently calculated for each of the two datasets share a common slope of 0.9 ‰ °C⁻¹ but differ in intercept by 7.2 ‰ (Fig. 6a). In the d-excess vs. δ¹⁸O relationship, under cold condition (i.e. for δ¹⁸O values lower than −45 ‰), the two sectors exhibit distinct compositions, with mean slightly higher d-excess values (∼ 2 ‰) observed for the Indian samples. Instead, in warmer coastal regions, d-excess shows higher variability (Fig. 6b). The comparison with bulk samples of EAIIST reveals that Section I – primarily influenced by Indian sector – exhibits isotopic values consistent with those of the corresponding divide of the Antarctic dataset. In contrast, Section II shows a pattern more aligned with Pacific sector.

Figure 7 (a) δ¹⁸O and (b) d-excess LMDZ6iso values vs. bulk observations. The black line is for linear regression.

We compare the observed and LMDZ6iso simulated bulk samples. For δ¹⁸O (Fig. 7a), the mean and standard deviation of the differences between observed and modelled δ¹⁸O values are −0.95 ± 2.24 ‰ (RMSE = 2.41 ‰), following well the different δ¹⁸O-T relationships between Indian and Pacific sector as observed in the Antarctic dataset (Appendix C, Fig. C1). For d-excess (Fig. 7b), the model overestimates the values for both sectors by approximately 5.5 ± 2.48 ‰ (RMSE = 6.07 ‰).

Figure 8 (a) δ¹⁸O and (b) d-excess composition of outward and return surface snow samples (dots) in comparison to the LMDZ6iso model outputs (lines). Vertical bars indicate spatial variability. (c) The days between the two journeys (grey) and the precipitation-free days prior the return sampling (black).

3.5 Temporal variability of surface isotopic composition

To assess the impact of post-depositional effects on the isotopic composition of the snow, we compare surface samples (representing a seasonal signal) collected during the outward and return ways of the EAIIST traverse across the plateau region (sampled between 10 and 50 d apart). The return isotopic composition is, on average, slightly higher for δ¹⁸O by approximately 3 ‰ (Fig. 8a), while d-excess is lower, by 5 ‰–10 ‰ (Fig. 8b). Limited snowfall for most of the sites over the plateau suggests that post-depositional processes could be responsible for the difference observed, allowing for prolonged interaction at the snow-atmosphere interface.

We evaluate the contribution of precipitation to the metamorphism of surface snow between the outward and return journey using the LMDZ6iso simulation outputs. The contribution of precipitation would lead to an increase of up to 2 ‰ for δ¹⁸O and a decrease of ∼ 1 ‰ for d-excess at sites north of 78° S. South of 78° S, the model suggests no change in surface snow isotopic composition due to precipitation. Using ERA5, we investigate the precipitation patterns between the outward and return samplings. In Section I, we identify 7–20 precipitation-free days and a total precipitation of 5–8 mm of fresh snow (corresponding to ∼ 25 % of the surface sample). In contrast, Section II experienced 10–20 dry days with no notable precipitation inputs (Fig. 8c).

3.6 Modelled sublimation effects on surface snow

In Sect. 3.5, we showed that precipitation input alone could not explain the change of isotopic composition of the surface samples between the outbound and return journeys. Here, we assess whether sublimation – commonly associated with an enrichment in δ¹⁸O along with a decrease of d-excess (Dietrich et al., 2023; Hughes et al., 2021) – can explain, at least in part, the variations observed. To this end, we compare the relative variations in d-excess and δ¹⁸O against the snow metamorphism conceptual model proposed by (Casado et al., 2021), which is able to link the relative change of d-excess and δ¹⁸O to the sublimated fluxes in the atmosphere (see Sect. 2.7). Overall, our observations support this framework: isotopic changes along Section II – during which no precipitation occur according to ERA5 – closely follow the modelled metamorphic modification of isotopic composition characterised by a slope of $\mathit{dxs}/{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ of −2  ‰ ‰⁻¹ (blue dots in Fig. 9) (Casado et al., 2021). For Section I, as shown in Sect. 3.5, precipitation accounts for a 25 % of the return sample weight, as a result, we expect a mix with the contibution of precipitation, characterised by a slope $\mathit{dxs}/{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ of −0.5  ‰ ‰⁻¹ (Casado et al., 2021). We observe that the surface snow in Section I evolved following a slope of −1.1  ‰ ‰⁻¹ (red dots in Fig. 9), which match the hypothesis of an input dominated by sublimation, with a contribution from the precipitation input.

Figure 9 d-excess vs. δ¹⁸O evolution of snow per day for Section I (empty red dots) and II (empty blue dots) centered on outward composition (black dot), in comparison with modeled $\mathit{dxs}/{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ sublimation slope (blue dashed line) and precipitation slope (black dashed line). The fluxes relative to sublimation are expressed in $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{d}}^{-\mathrm{1}}$. The errorbars represent the uncertainty from spatial variability.

Figure 10 Difference between return-snow observations and modelled δ¹⁸O and d-excess values, for Section I (red) and Section II (blue). Modelled values are calculated considering either precipitation only (crosses) or both precipitation and sublimation effects (circles). Error bars represent the uncertainty of the LMDZ6iso model in simulating the isotopic precipitation at Concordia (Dutrievoz et al., 2025). The x-axis indicates the cumulative precipitation between the outward and return samplings.

We use Section II data, driven solely by sublimation, to estimate the flux of water transfered from the snow to the atmosphere (eq. Sect. 2.7). The estimated sublimation fluxes range from 0.04 to 0.09 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{d}}^{-\mathrm{1}}$ (Fig. 9), falling towards the lower end of previous estimates (0.05–0.35 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{d}}^{-\mathrm{1}}$; Ollivier et al., 2025). It is important to note that the values reported by Ollivier et al. (2025) represent daily estimates and therefore reflect substantial day-to-day variability driven by changes in insolation, cloud cover, temperature, and wind speed, as well as processes such as nighttime condensation. In contrast, our estimates represent averages over periods of up to 50 d. Therefore, they integrate over both days with intense sublimation and days with limited or no sublimation (e.g., cloudy or snowy conditions), resulting in a narrower overall range.

We use this sublimation impact for modeling the final composition of the snow along the whole traverse, using outward snow samples as the initial condition. This calculation accounts for the number of precipitation-free days between samplings (from ERA5), and incorporates precipitation input from the LMDZ6iso model. The final δ¹⁸O and d-excess modeled composition is estimated as:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\text{Final\hspace{0.17em}\hspace{0.17em}Snow}}_{\text{model}}=& {\text{Initial\hspace{0.17em}\hspace{0.17em}Snow}}_{\text{obs.}}\\ & +{\text{Precipitation\hspace{0.17em}\hspace{0.17em}Input}}_{\text{LMDZ6iso}}\\ & +(\text{Sublimation}\times \mathrm{N}{\mathit{°}}_{\text{dry\hspace{0.17em}\hspace{0.17em}days}})\end{array}\end{array}$$

The difference between observed versus modelled final values are shown in Fig. 10. Red symbols represent Section I sampling sites, characterized by high precipitation (> 1 $\mathrm{mm}\mathrm{w}.\mathrm{e}.$) between outbound and return samplings, for which freshly precipitation represents ∼ 25 % of the sampled snow. Blue symbols represent Section II sites with low precipitation (< 1 $\mathrm{mm}\mathrm{w}.\mathrm{e}.$), characterised by negligible precipitation. The modelled values are presented considering either only the precipitation input or both precipitation and sublimation effects. Including sublimation in the computation reduces the discrepancy for all Section II sites and for the majority of Section I sites. The mean absolute difference decreases from 1.9 to 1.3 ‰ for δ¹⁸O, and from 6.6 to 2.9 ‰ for d-excess.

3.7 Discussion

3.8 Air mass origins and δ¹⁸O-temperature slopes

The spatial δ¹⁸O-T relationship has traditionally been used to calibrate isotopic paleothermometer (Dansgaard, 1964; Lorius et al., 1969). However, several studies have shown that it is generally less accurate than the δ¹⁸O-T temporal slope derived at a single site (Sime et al., 2008; Casado et al., 2017). The EAIIST traverse, which crosses two sectors of East Antarctica influenced by distinct air mass origins yet characterised by relatively similar climatic conditions, provides valuable insights into the role of moisture sources in shaping both spatial and temporal slopes. These observations also highlight the potential to model such relationships in regions where long-term temporal records are not available.

Figure 11 δ¹⁸O-T relationship for monthly samples simulated by LMDZ6iso at Dome C, AGO5, PALEO and MD sites. The color scale indicates the dominant longitudinal origin of air masses associated with each monthly precipitation event. Significant linear relationship (p value < 0.05) is indicated by asterisks.

The spatial analysis of EAIIST Section I showed a significant δ¹⁸O-T relationship over the Indian sector (p < 0.05; Table 1), while no significant correlation was detected for Section II. It is important to note that Section II has limited spatial coverage, as the traverse did not extend to the Pacific coast, and the temperature range along this 300 km transect is correspondingly narrow. Therefore, the absence of a detectable correlation cannot be interpreted as evidence that a δ¹⁸O–T relationship is lacking in this sector, as confirmed by analysis of the broader Antarctic dataset (Fig. 6a), which shows statistically significant correlations in both sectors (p < 0.05). The spatial slopes for the two basins are similar (0.9 ‰ °C⁻¹), reflecting comparable thermal history of air mass from the coasts to the high-elevation regions in the interior of the continent (Helsen et al., 2007), while differences in the intercepts are attributed to contrasting vapor isotopic compositions associated with distinct source regions and longer transport pathways for Pacific air masses.

Results from LMDZ6iso model reproduce the inter-sector differences in δ¹⁸O-T spatial slopes for EAIIST (Figs. 7a and C1a), and for the broader Antarctic dataset (Fig. C2). Although a comprehensive quantitative evaluation of the LMDZ6iso model is beyond the scope of this study, the qualitative comparison indicates that the model successfully captures the observed spatial variability associated with differing air mass influences. This supports the hypothesis that LMDZ6iso is also able to simulate temporal δ¹⁸O-T slopes. However, given the limited spatial coverage of Section II, future field campaigns should be targeted toward the Pacific sector on the East plateau to further investigate and validate this hypothesis.

Monthly δ¹⁸O-T relationships, derived from modelled precipitation at single sites (i.e. Dome C, AGO5, PALEO and MD; Fig. 11), account for varying moisture sources and yield temporal slopes of 0.4–0.5 ‰ °C⁻¹ (R² = 0.9). At Dome C, this prediction is consistent with observed monthly precipitation collected at Concordia Station, which exhibits a temporal slope of ∼ 0.5 ‰ °C⁻¹ (Dreossi et al., 2024).

We note that here the term “snow-precipitation samples” refers to both model outputs and precipitation samples collected immediately after deposition where snow was captured on a wooden platform positioned 1 m above the surface to minimize post-depositional alterations. In contrast, the term “snow-surface sample” refers to material collected directly from the surface during the traverse, representing an integrated signal of multiple precipitation events and influenced by wind redistribution and post-depositional metamorphism.

3.9 Modelled d-excess vs. δ¹⁸O relationship

The d-excess vs. δ¹⁸O relationship (Fig. 6b) shows high d-excess variability in warmer coastal areas (< 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$), reflecting sensitivity to local moisture-source conditions such as temperature, relative humidity in evaporative region, and sea ice extent. Inland (> 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$), at the end of the distillation pathways, d-excess becomes more sensitive to condensation temperatures, since it is regulated by independent supersaturation functions (Uemura et al., 2012; Stenni et al., 2016; Touzeau et al., 2016; Landais et al., 2017). This pattern is evident in the Antarctic database and is weaky detectable in EAIIST for very cold sites (δ¹⁸O < −50 ‰), where Indian-sector samples show slightly higher d-excess than Pacific-sector samples at similar δ¹⁸O values.

LMDZ6iso model does not accurately predict the d-excess values, showing a RMSE equal to 6.07 ‰ (Figs. 7b and C1b). The mismatch primary arises from the model formulation of the isotopic supersaturation parameter (λ) in the AGCM. The chosen value of λ equal to 0.004 K⁻¹ is selected to achive a optimal compromise to correctly simulate both δ¹⁸O and d-excess in surface snow, but leads to a systematic overestimation of d-excess (Dutrievoz et al., 2025). In addition, the model do not include post-depositional processes. Sublimation lowers d-excess in surface snow by modifying isotopic composition after deposition (Landais et al., 2017; Wahl et al., 2022). This effect is particularly evident in the observations at Windcrust site (WC; Fig. 5), where low d-excess values (∼ 4 ‰) reflect strong wind activity, substantial erosion of accumulation, and enhanced snow–atmosphere exchange. The influence of wind-snow redistribution and sublimation on the isotopic signal is discussed in detail in Sect. 4.3.

3.10 Impact of post-depositional processes

On the Antarctic plateau, post-depositional effects mainly include wind-driven snow redistribution, and sublimation.

Snow transport by wind is a relatively local process, typically mixing snow from the surrounding areas and one of the main contributions to stratigraphic noise (Hirsch et al., 2023). Studies suggest that the snow shuffled by wind redistribution can reach distances up to ∼ 100 km (Scarchilli et al., 2010; Frezzotti et al., 2007). Such mixing generates stratigraphic effect between precipitation events, leading to δ¹⁸O variability of up to 4.4 ‰ within the uppermost 6 cm at Kohnen Station in the plateau interior (Münch et al., 2016). In our dataset, this effect is reduced because we mixed surface snow collected over an extended area at each sampling site and by integrating 1 m of snow depth for bulk samples.

To quantify sublimation, mainly impacting the snow metamorphism in summer, we compared surface snow samples collected during the outbound and return ways of the traverse, where low precipitation and strong snow drift favour prolonged summer atmosphere-snow exposure. Assuming 20 precipitation-free days per summer season, we estimated a sublimation-induced variation corresponding to +2 ‰ yr⁻¹ in δ¹⁸O and −8 % yr⁻¹ in d-excess (see Sect. 3.6).

Including this sublimation effect in the modelled isotopic predictions substantially improved the agreement with observations, reducing the discrepancy for d-excess compared to simulations considering precipitation input alone (Fig. 10).

This improvement is particularly evident at sites where cumulative precipitation between samplings is negligible (< 1 $\mathrm{mm}\mathrm{w}.\mathrm{e}.$). The remaining differences and variability between observed and modelled values can be partially attributed to uncertainties in ERA5 precipitation and LMDZ6iso model. Previous studies have shown ERA5 overestimates precipitation over the East Antarctic Plateau, with biases reaching up to 50 % relative to satellite-based measurements (Roussel et al., 2020). As a result, our modelling likely represents the maximum contribution of precipitation, implying that the metamorphism would be even greater if actual precipitation were lower. We emphasize that the aim of this study is not to quantify the ability of the combined ERA5-LMDZ6iso in reproducing the absolute isotopic values, but rather to evaluate whether accounting sublimation improved the qualitative representation of surface snow isotopic composition compared to precipitation-only scenarios. These results further reinforced the key role of post-depositional processes in shaping the isotopic composition of surface snow.

4 Conclusions

In this study, we presented new observational and modelling evidences on the processes shaping the isotopic composition of surface snow across the East Antarctic Plateau. Together, these results provide key insights for the future interpretation of ice cores collected in low-accumulation area:

• i.

  Observational snow samples defined the range of isotopic variability associated with distinct moisture-source sectors, an essential consideration in regions where multiple air-mass trajectories coexist and may shift over time.

• ii.

  The comparison between observations and LMDZ6iso outputs confirmed the model's ability to reproduce both spatial and temporal δ¹⁸O-T slopes. The temporal slope, in particular, is crucial for calibrating the paleothermometer in areas lacking long-term records.

• iii.

  Sublimation has a great influence on the isotopic signal, particularly decreasing d-excess in summer. This demonstrates that achieving a robust qualitative assessment of isotopic variability requires the explicit inclusion of isotopic fractionation during sublimation and condensation processes in General Circulation Models (GCMs).

• iv.

  Finally, the finding that relatively close drilling sites (e.g. Dome C and PALEO, separated by only few hundred kilometres) may record precipitation from distinct moisture pathways suggests a promising strategy to enhance temperature reconstructions combining ice core records that capture common large-scale climate variability through different precipitation events.

Appendix A: Londitudinal distribution of air mass orising on the East Antarctic Plateau sampling sites

Figure A1 Longitudinal distribution back-trajectories over the 2009–2019 period for four representative sites: Dome C, AGO5, PALEO, and MD. The median longitudinal origin of air masses for the 33 sampling sites across the Antarctic Plateau is indicated by black dots.

Appendix B: Spatial distribution of surface and bulk samples

Figure B1 Spatial distribution of δ¹⁸O and d-excess values of surface (empty markers) and bulk (filled markers) samples for Section I (red – north of 78° S) and Section II (blue – south of 78° S). Temperature, elevation and distance from DDU are reported.

Appendix C: LMDZ6iso predictions

Figure C1 (a) δ¹⁸O-T relationship and (b) d-excess vs. δ¹⁸O for LMDZ6iso bulk samples of the traverse compared to observations from Antarctic database, separated into Pacific and Indian sectors.

Figure C2 δ¹⁸O-T spatial relationship predicted by LMDZ6iso for the Antarctic dataset, separated into Pacific and Indian sectors.