Water Isotopic Signature of Surface Snow Metamorphism in Antarctica

Water isotope ratios of ice cores are a key source of information on past temperatures. Through fractionation within the hydrological cycle, temperature is imprinted in the water isotopic composition of snowfalls. However, this signal of climatic interest is modified after deposition when snow remains at the surface exposed to the atmosphere. Comparing time series of surface snow isotopic composition at Dome C with satellite observations of surface snow metamorphism, we found that long summer periods without precipitation favor surface snow metamorphism altering the surface snow isotopic composition. Using excess parameters (combining D,17O, and 18O fractions) allow the identification of this alteration caused by sublimation and condensation of surface hoar. The combined measurement of all three isotopic compositions could help identifying ice core sections influenced by snow metamorphism in sites with very low snow accumulation.

understood post-deposition processes that cannot be isolated and removed (Petit et al., 1982;Waddington et al., 2002), and which may affect the temperature reconstructed from water isotopic composition along the ice core.
The exchange between the firn and the atmosphere is likely to override the initial climate signal in the precipitating snow, and most post-deposition processes lead to such exchanges. First, the cumulated contributions of sublimation and surface condensation can exceed 10% of the annual surface mass balance  which can affect the snow isotopic composition signal (Casado et al., 2018;Madsen et al., 2019;Steen-Larsen et al., 2014). Second, surface snow, exposed to solar radiation and the atmosphere, is subjected to strong metamorphism, which is due to the microscale transport of vapor between the snow grains and interstitial air spaces (Wiscombe & Warren, 1980). Third, wind pumping can lead to vapor advection in and out of the snowpack at depth (i.e., the firn) (Calonne et al., 2014). All these processes have been reported to cause isotopic exchanges between the surface snow, the atmosphere above (Sokratov & Golubev, 2009;Steen-Larsen et al., 2014;Town et al., 2008) and the interstitial air below the surface (Ebner et al., 2017;Johnsen, 1977;Neumann & Waddington, 2004;Petit et al., 1982). Overall, the isotopic signal stored in ice cores can be seen as the combination of the following contributors (Casado et al., 2018(Casado et al., , 2020Laepple et al., 2018): (a) temperature signal, (b) noise linked with precipitation intermittency, (c) signal from post-deposition processes (including snow metamorphism), and (d) isotopic diffusion within the firn. The effects of post-deposition processes are mainly documented on time scales from hours to days (Casado et al., 2018;Hughes et al., 2021;Ritter et al., 2016;Steen-Larsen et al., 2014) and then at interannual time scales, and on the East Antarctic Plateau, the climatic signal is expected to account for only 10%-50% of the variance in δ 18 O (Casado et al., 2020;Laepple et al., 2018;Münch & Laepple, 2018).
Here, we use the combination of water isotopes (δD, δ 18 O, and δ 17 O) in the second order parameters d-excess and 17 O-excess to explore the isotopic signature of post-deposition processes associated with metamorphism in surface snow. Because of the different sensitivities of δD, δ 18 O and δ 17 O to equilibrium and kinetic fractionation, d-excess and 17 O-excess are indeed expected to be strongly influenced by water vapor diffusion associated with surface metamorphism. We explore the link between the surface snow isotopic composition and grain size, the latter being a measure for the degree of exposure to metamorphic processes at the surface (Ebner et al., 2017): During metamorphism, sublimation and re-condensation on slightly colder grains lead to a coarsening of the snow grains. This process is driven by a surface temperature gradient exceeding 5 K m −1 affecting the first few centimeters of snow in Antarctica (Picard et al., 2012). Such large surface temperature gradients are expected to create moisture fluxes (Gallet et al., 2014) which affect the water isotopic repartition in surface snow (Casado, 2018;Landais et al., 2017).

Surface Snow Isotopic Composition
The sampling and measurement of surface snow isotopic compositions is detailed in (Casado et al., 2018). Bi-weekly duplicate snow samples were collected between December 2013 and April 2018 at Dome C, Antarctica, at random locations within a 100 × 100 m area. For δ 18 O and d-excess, 406 samples were measured by infrared spectroscopy while 17 O-excess was measured by mass spectroscopy with a limited throughput, only 72 were analyzed (Figure 1). The 95% confidence intervals are calculated by computing a standard deviation of the difference between replicates on 72 data points (datapoints where all three isotopes were available, for δ 18 O and d-excess, we obtain similar results using the ensemble of 406 samples).

Numerical Evaluation of the Precipitation Input
We simulate the surface snow isotopic composition in the case of precipitation input only (Casado et al., 2018). To do so, we built an evolving profile of δ 18 O and isolate for each day the isotopic composition of the top 1.5 cm. The δ 18 O profile is built by adding at the surface, for each precipitation event, a new layer of snow for which its isotopic composition is obtained by converting the temperature: The factors are obtained from comparing precipitation isotopic composition and temperature at Dome C (Stenni et al., 2017;Touzeau et al., 2018). The thickness of the new layer is given by the amount of precipitation of the event and the snow density (set constant at 180 kg m-3 [Leduc-Leballeur et al., 2017]). The same exercise is done iteratively for each precipitation event using the temperature and precipitation fields from ERA-interim.
To include surface snow d-excess in this model, we use the slope between d-excess and δ 18 O of 1.4 ‰/‰ found in precipitation isotopic composition . Sensitivity tests were realised to see if a lower slope between d-excess versus δ 18 O affects the results (Supporting Information S3).

Numerical Evaluation of the Metamorphic Input
The metamorphic input includes two different processes: Sublimation, occurring at dome C mainly during the warmest months of the year, and condensation, occurring when the surface becomes colder than the atmosphere or the sub-surface.
In our model, during periods dominated by sublimation, the isotopic composition of surface snow remaining after sublimation i snow R under the local closure assumption and in steady state (Dongmann et al., 1974;Merlivat & Jouzel, 1979) is: i eq sub corresponds to the equilibrium fractionation coefficient associated with sublimation with the snow; RH is the air relative humidity (as defined by [Merlivat & Jouzel, 1979], see Supporting Information S4a); and i a R is the isotopic composition of the atmospheric vapor. The exponent of k is the classical exponent used in the Craig and Gordon equation, which accounts for the roughness of the boundary layer turbulence (Merlivat & Jouzel, 1979). To apply Equation 2 to our problem, we assumed a relative humidity of 80% with respect to the liquid saturated vapor pressure, which corresponds to a partial pressure of 65 mbar, a typical value for summer at Dome C (Andreas et al., 2002;Genthon et al., 2013;Gettelman et al., 2006). We explored the space of parameters (different sets of fractionation coefficients and terms in Equation 2), and observed that the roughness parameter k = 0.4 fits best the observations (See Supporting Information S4a).
The isotopic composition of the condensate is obtained as in the Crocus-iso model (Touzeau et al., 2018): where   i eq cond and   i eq sub are the equilibrium fractionation coefficients associated with condensation at the surface and sublimation at the sub-surface where the moisture originates from, respectively, and  [Casado et al., 2018]), δ 18 O (lighter green: observations, darker green: precipitation modeling approach [Casado et al., 2018]), grain index (black) (Picard et al., 2012), temperature from ERA-interim reanalysis (red), and precipitation amounts from ERA-interim reanalysis (gray bars). For the isotopic compositions the shaded area indicates uncertainty range (see Methods). The vertical shading indicates the periods when intense metamorphism is detected using the grain index, a red bar indicates a period during which metamorphism is associated with sublimation (Phase 1 below) while a blue bar indicates a period associated with condensation (Phase 2 below) (b) schematic representation of the processes involved during the sublimation phase; (c and d): Evolution of the isotopic composition linked to each elementary process in the space δD versus δ 18 O and ln (δ 17 O + 1) versus ln (δ 18 O + 1), respectively; while panels (b-d) show Phase 1, (e-g) show the equivalent figures for Phase 2 (sublimation in red, condensation in blue).
By including the input of moisture from the sub-surface, our model could thus reconcile the large difference between the observed mass gain at the snow surface and the small amount of moisture being deposited from the atmosphere (Supporting Information S1). During this period, there is very little vapor in the atmosphere. Nevertheless, we predict that while the atmospheric contribution in terms of the total number of molecules is insignificant, the atmosphere still exchanges isotopes with the surface snow and the interstitial vapor which creates conditions for the condensation at the snow surface that are out of equilibrium. Condensation of surface hoar is a process which mainly involves vapor diffusion, and for which, we expect kinetic fractionation to play a major role for the isotopic content of the surface snow.

Time Series of Surface Snow Isotopic Composition
The classical interpretation of the influence of temperature on δ 18 O considers that surface snow δ 18 O is higher during warm periods and lower during cold periods, which is in overall agreement with our measurements obtained at Dome C ( Figure 1a ). This is also the case for surface snow at Dome C (see Figure 1a) and Supporting Information S1). However, the negative correlation observed in surface snow isotopic compositions (r = −0.51, p < 0.05) is weaker than expected compared to spatial transect across Antarctica (Masson-Delmotte et al., 2008) and the amplitude of the relative variations in d-excess and δ 18 O show strong variations: In 2013-2014, when large amounts of precipitation were deposited (twice the long-term average), we observe small variations in d-excess (a decrease of roughly 10‰) coincident with large variations in δ 18 O (an increase of 18‰, see Table S1); whereas in 2014-2015, when barely any precipitation was deposited, we observe large variations in d-excess (a decrease by up to 20‰) coincident with small variations in δ 18 O (an increase of only 6‰).
Similarly, we observe that the surface snow 17 O-excess is negatively correlated to δ 18 O during summer 2014-2015 (r = −0.26, p < 0.05). This is opposite to what is expected from fractionation coefficients, observations of precipitation isotopic composition, and from long-term ice core records ( Figure S2), but similar to observations taken at another site with very low accumulation on the East Antarctic Plateau: The Vostok station in a three-meter snowpit (Table S1) (Landais et al., 2012;. The main 17 O-excess variation is larger (28 ppm between the end of January and mid-February) during the 2014-2015 summer than during the 2013-2014 summer (15 ppm), which is similar to the signal observed for d-excess (Figure 1a).

Signal Explained by the Precipitation Input
We evaluate which part of the surface snow isotopic composition signal can be explained by the precipitation input, as classically interpreted. To do so, we apply the modeling approach solely based on the precipitation input (See Section 2.2) to showcase the other contributions to the surface snow isotopic composition (Casado et al., 2018). For the summer of 2013-2014, more than 70% of the observed changes of surface snow δ 18 O (p < 0.05, N = 36) can be explained by the dominant input of precipitation (dark green line in Figure 1). By contrast, for the summer 2014-2015, the surface snow δ 18 O variance explained by the precipitation input is r 2 = 0.02 (p = 0.4, N = 36). This result is line with studies using vertical profiles of firn which found that more than 50%-90% of the variability at the seasonal scales was created by surface processes (Casado et al., 2020;Laepple et al., 2018;Münch & Laepple, 2018). Not only is the amplitude of the simulated changes at odds with the observations, we show a two-month shift between calculated and measured variations (Figure 1a).
The strong differences between water isotopic composition predicted by the precipitation modeling approach for each phase of summer 2014-2015, as well as the strong difference of the relationship between δ 18 O and d-excess in 2013-2014 and 2014-2015 argue for the need for contributions other than solely precipitation to understand surface snow isotopic composition. We evaluate the relationship between surface snow δ 18 O and temperature for these two consecutive years: 2013-2014 when precipitation seemed to be the main driver of the isotopic signal and 2014-2015 when metamorphism could play a significant role. As the surface snow integrates several events of precipitation, it is not possible to calculate the δ 18 O versus temperature relationship with a simple linear regression because of the hysteresis behavior of the two variables (Casado et al., 2018). Here, the relative changes of surface snow δ 18 O versus temperature in summer 2013-2014 are 0.49‰/°C (Casado et al., 2018), very close to the value obtained for the precipitation δ 18 O versus temperature (0.46‰/°C ). Applying the same method to evaluate the relationship in 2014-2015 than in (Casado et al., 2018), we obtain a relationship of 0.24‰/°C, much lower than the signature of the precipitation input.

Additional Contribution of Snow Metamorphism
We compare the variations of surface snow isotopic composition with the ones of snow physical properties by using the grain index, a satellite product that evaluates the average grain size in the surface top 7 cm (Picard et al., 2012). A much smaller grain index signal was observed in the summer 2013-2014 (from 01/12/2013 to 15/01/2014, increase of 0.07) compared to summer 2014-2015 (from 01/12/2014 to 15/01/2015, increase of 0.17). The summer associated with low grain index increase (2013-2014) can be explained by large precipitations events, which bring a significant amount of small grains of snow (Picard et al., 2012). In contrast, it has been shown that during the summer 2014-2015, the snow remained exposed at the surface for a long time, facilitating a higher degree of metamorphism, which led to a large increase of the grain size (Leduc-Leballeur et al., 2017).
In both summers, the evolution of δ 18 O shows an increasing and decreasing phase (Figure 1a). For the summer of 2013-2014, the δ 18 O increasing phase (Phase 1) lasts from November 25, 2013 to January 27, 2014 and is associated with relatively small decreases in d-excess and 17 O-excess compared to the δ 18 O increase of 18‰ (slopes of −0.22‰/‰ and −0.6 ppm/‰ (non-significant), respectively, see Table S1). This is followed by the Phase 2 with a progressive and steady decrease in δ 18 O until April 15, 2014. In the summer of 2014-2015, Phase 1 spans a similar period (December 4, 2014 to January 29, 2015) but features only a 6 ‰ increase in surface δ 18 O (Figure 1). In parallel, d-excess and 17 O-excess decreased by 12 ‰ and 12 ppm, respectively, associated with a negative correlation (slopes of −2.0 ‰/‰ and −2 ppm/‰, respectively, see Table S1). Phase 2 lasted from January 30, 2015 to March 20, 2015 with δ 18 O slightly decreasing (by 2.5 ‰) while both d-excess and 17 O-excess exhibit large increases of 12‰ and 23 ppm, respectively. This period is characterized by a significant decrease in air temperature and ends on March 20, 2015 with an abrupt drop in the grain index, caused by an intense wind event eroding a 2 cm surface layer (Leduc-Leballeur et al., 2017).

Discussion
The concomitant changes of surface snow isotopic composition and grain index suggest that metamorphism could be the source of the observed anomalous signal in summer 2014-2015. We therefore developed simple models (see Methods) to explain the anomalous variations in δ 18 O, d-excess, and 17 O-excess observed during the summer of 2014-2015 when no significant precipitation events occurred.

Phase 1: Sublimation
Our data clearly show that we have isotopic fractionation in this period which is otherwise characterized by sublimation , in opposition to the classically assumed process of uniform sublimation of ice layers on perfectly spherical snow grains in line with (Madsen et al., 2019) and (Wahl et al., 2021). This is likely a result of the porosity of the snow grains, for which it was shown that sublimation would occur with isotopic fractionation (Ebner et al., 2017;Hughes et al., 2021;Touzeau et al., 2018).
Our steady-state model (see Methods) is inspired by the Craig and Gordon approach (Craig & Gordon, 1965), the isotopic composition of the formed water vapor is calculated taking into account the fractionation at (a) the phase transition (here sublimation), (b) diffusion through the boundary layer, and (c) exchange CASADO ET AL.

10.1029/2021GL093382
with the free atmosphere (Figure 1b)). We use a schematic representation of sublimation in our model from (Equation 2) and separate the three fractionation processes to determine their respective influence on d-excess and 17 O-excess (Figure 1c and 1d). First, equilibrium fractionation above the surface snow (term 1 in Equation 2) occurs with slopes slightly higher than the meteoric water lines (see arrow 1 in Figures 1c  and 1d). This process leads to a large increase in δ 18 O and small increases in d-excess and 17 O-excess. The second process is the kinetic fractionation associated with diffusion with a slope well below the meteoric water line (term 2 in Equation 2). This process leads to a further, yet moderate, increase of δ 18 O and a strong decrease of both d-excess and 17 O-excess during Phase 1 (arrow 2 in Figures 1c and 1d). The third process is the mixing with the atmospheric vapor (term 3 in Equation 2, represented by arrow 3 in Figures 1c and 1d) which leads to a decrease in δ 18 O due to dilution with the δ 18 O-depleted local atmospheric vapor.
In contrast to the signals obtained when considering the precipitation modeling approach, our model taking into account sublimation is able to quantitatively reproduce the negative correlation between d-excess (and 17 O-excess) (Table S1). We found that the isotopic composition of the free atmosphere ( i a R ) is a crucial parameter that influences the changes of 17 O-excess of the surface snow isotopic composition during the sublimation (see sensitivity tests in the Supporting Information S4a). The roughness of the turbulence in the atmosphere is also a key parameter (Table S3), in the case of a smooth regime (k > 2/3 instead of 0.4), diffusion would induce a larger isotopic fractionation, leading to increased d-excess and 17 O-excess signals. Weakly turbulent boundary layer could lead to such signals, but do not seem to have occurred during the period of our observations.

Phase 2: Condensation of Surface Hoar
The formation of hoar observed during Phase 2 of the summer 2014-2015(Leduc-Leballeur et al., 2017 was likely due to the rapidly decreasing temperature of the surface linked to decreased radiative input. When this happened, the snow surface became colder than both the subsurface layer and the atmosphere ( Figure S2). During this period, only one precipitation event accounting for 0.5 kg m −2 occurred, while the amount of mass gained at the surface was 3 ± 2 kg m −2 (see Supplementary Material). This large difference cannot be explained by condensation from the atmosphere because moisture fluxes from the atmosphere toward the snow were estimated not to exceed 0.05 kg m −2 (see Supporting Information S1). Our model calculates the isotopic fractionation due to condensation at the surface of snow sublimated at the subsurface (See Methods).
The condensation process can be divided into 3 steps (Figures 1f and 1g). First, we consider the water vapor formed at isotopic equilibrium with the snow at the sub-surface (term 1 in Equation 3). Then, as the sub-surface temperature is higher than the surface temperature, lighter isotopes in the interstitial vapor are diffusing preferentially toward the colder surface snow (term 2 in Equation 3). This diffusion is followed by condensation of supersaturated water vapor with another step of equilibrium fractionation (term 3 in Equation 3). Overall, these three steps lead to a decrease in δ 18 O and an increase in both d-excess and 17 O-excess since diffusion is associated with slopes lower than the meteoric water lines in δD versus δ 18 O and δ' 17 O versus δ' 18 O spaces. The changes of δ 18 O, d-excess, and 17 O-excess are in qualitative agreement with the evolution of the isotopic composition observed in the field during Phase 2 of the summer of 2014-2015 (Table S1). Based on sensitivity tests performed for our different parameters, the low turbulence solution best fits the observations, which is realistic considering diffusion inside the porosity of the snow (Supporting Information S4b). The main uncertainty of our calculation comes from the isotopic composition of the sub-surface snow, which remains as a free parameter. Without observation to constrain the sub-surface snow isotopic composition, only qualitative evaluations can be made (see Tables S9 and S10 compared to  Table S8).

Generalization of the Metamorphic Isotopic Signature
We identified that three of the contributors to the surface mass balance (SMB) namely precipitation, sublimation, and condensation, contribute to the surface snow isotopic content with different isotopic signatures (Table S1 [Hughes et al., 2021]). We were able to characterize these signatures making use of two particular summers, one associated with a dominating precipitation input to the SMB (summer 2013-2014, during which Δd-exc/Δδ 18 O = −0.5 ‰/‰), and one associated with a dominating metamorphism (summer 2014-2015, Δd-exc/Δδ 18 O between −2.5 and −5.5 ‰/‰), and to show that the relative changes observed can be explained by a simple modeling approach of the different fractionation processes.
We extend the results of the study of the surface snow isotopic composition signature from December 2013 to April 2018 to evaluate the relevance of metamorphism influence on ice core isotopic composition (Figure 2a). Over this extended period, the precipitation input can only explain 48% of the total variance of the surface snow isotopic composition variations (p < 0.05, N = 373). During the periods dominated by precipitation, we use two different parametrisations: A d-excess versus δ 18 O slope of −1.4 ‰/‰, similarly to the one used in the precipitation modeling approach, and a d-excess versus δ 18 O slope of −0.5 ‰/‰, as found during the summer 2013-2014 when no impact of metamorphism has been identified (See Supporting Information S5). We generalize the finding that the surface snow isotopic composition cannot be predicted by the sole precipitation input (Figure 2a), light and dark green and blue lines for δ 18 O and d-excess, respectively): Differences between the predicted and observed δ 18 O reach 8, 2.5 and 7.5‰ in summer 2014-2015, 2015-2016 and 2017-2018, respectively, while for d-excess, the differences over the same summers are up to 16, 9, and 7‰. Over the extended period (2013)(2014)(2015)(2016)(2017)(2018), the dominant mode of variability for the surface snow isotopic composition (PC1 calculated from the δ 18 O and d-excess time series, Figure 2a), orange line, associated with a slope of −1.6 ‰/‰), explaining 78% of the variance, suggest that both the precipitation and the metamorphic inputs play a significant role in the surface snow isotopic composition signal at Dome C ( Figure 2b). Indeed, the PC1 captures both the strong d-excess negative peak of summer 2014-2015 linked with metamorphism and the δ 18 O signal of 2013-2014 linked to precipitation, and its isotopic signature ranges between the one identified for the precipitation (between −0.5 and −1.4 ‰/‰) and the one identified for the metamorphism (between −2.0 and −4.9 ‰/‰).

Metamorphic Isotopic Signature Archived in Ice Core Records
As for the surface snow, the precipitation input does not explain a large part of the variability observed in the Vostok snowpit Figure 2c. Interestingly, the largest discrepancies between the snowpit and the precipitation modeling approach occur simultaneously with the larger values of d-excess and 17 O-excess, as CASADO ET AL.  well with the strong anti-correlation between d-excess ( 17 O-excess) and δ 18 O. These strong d-excess and 17 O-excess signals suggest that post-deposition processes are needed to be taken into account to explain the isotopic signal in ice cores, beyond the climatic signal imprinted in precipitation isotopic composition, or even precipitation intermittency (Casado et al., 2020), and justify using the same PCA diagnostic on the snowpit series. We performed the same principal component analysis for d-excess and δ 18 O than for the surface snow on snow pits from Dome C ( Figure S4A and S4B, ) and surrounding sites (Point S2, 160 km away from Dome C, Figure S4C and S4D, , and Vostok Figure 2c and 2d ).
The principal component analysis also highlights that the dominant mode of variability in the Vostok snowpit (75% of the variance explained) is associated with an isotopic signature of Δd-exc/Δδ 18 O = −3.2 ‰/‰, suggesting that a significant part of the signal is associated with a metamorphic input. Similar results are found with an ice core drilled at the point S2 ( Figures S4C and S4D, slope of −1.3 ‰/‰, 91% of explained variance) and to a lesser extent for an undated snow pit from Dome C ( Figures S4A and S4B, slope of −1.3 ‰/‰, 62% of explained variance). The PCA reveals that the metamorphism influence on δ 18 O and d-excess at the surface of the firn is preserved in the archived snow.

Conclusions
We showed that d-excess and 17 O-excess in snow at Dome C, Antarctica, are strongly influenced by summer metamorphism because the precipitation amount is low and the existing snow remains exposed to the atmosphere for a long time. By comparing the variations in δ 18 O, d-excess and 17 O-excess, it could be possible to best document the processes that drive SMB changes due to precipitation, sublimation toward the atmosphere, or condensation from underlying firn.
These results have consequences for climate reconstructions based on ice core records, because metamorphism can lead to additional isotopic signal that may remain imprinted on the ice isotopic composition. The 28 ppm increase in 17 O-excess of surface snow at Dome C during the condensation period (Phase 2 of the summer of 2014-2015) is equivalent to three times the signal observed during the last deglaciation (Winkler et al., 2012). As the intensity of both the sublimation and condensation processes shows a considerable seasonal and interannual variability (dependent on radiation, and cloud cover), as well as a variability at longer time scales (e.g., the glacial-interglacial), these processes could introduce biases in the interpretation of ice cores. Indeed, the link between δ 18 O and temperature when the precipitation input dominates (0.46 ‰/°C) is very different than the one associated with prevailing metamorphic influence (0.24 ‰/°C).
Another important application of data on post-deposition impacts on snow isotopic composition is to examine past variations in firn processes. Since our model can quantitatively reproduce the surface snow isotopic composition, measurements of 17 O-excess and d-excess could be used to quantify the importance of surface hoar formation, thereby complementing density or grain index measurements to determine local metamorphism. We suggest that isotopes can be used as a tracer of surface metamorphism to yield additional insights into the physics. Our simple model provides a new set of equations to simulate fractionation associated with sublimation and condensation which can be used in isotope enabled GCM and snow pack models to help quantifying metamorphism in the past from ice core records.