Quantifying the Stable Water Isotopologue Exchange Between the Snow Surface and Lower Atmosphere by Direct Flux Measurements

Surface processes in high latitudes play an important role in global climate and thus understanding the physics of these systems is critical for improving climate projections. The characterization of the stable water isotopologue flux between the surface and the atmosphere offers the potential to constrain parameterizations of these physical surface exchange processes in numerical models. In addition, observations of isotopologue surface fluxes allow the evaluation of surface fluxes as a process influencing the formation of the climate signal retrieved from ice core isotopologue records. Here, we present a record of stable water isotopologue surface fluxes measured in‐situ in the accumulation zone of the Greenland Ice Sheet at the East Greenland Ice Core Project site. We measured isotopologue fluxes above the snow surface directly by combining high‐frequency eddy covariance measurements with low‐frequency isotopologue measurements from a cavity ring‐down spectrometer (CRDS). We developed a method to correct for the high‐frequency loss of the CRDS by combining humidity measurements from both the CRDS and eddy covariance instruments. Using this approach our measurements provide the first direct observations of water isotopologue fluxes in polar areas. We observed a clear diurnal cycle in the fluxes of the different water isotopologues. The isotopic composition of the sublimation and deposition flux showed to be dependent on the snow and vapor isotopic composition, respectively. To a first order, the isotopic composition of the sublimation flux could be derived assuming equilibrium fractionation during sublimation.

Greenland Ice Sheet (GrIS) helps to maintain the total mass of the ice sheet (Box & Steffen, 2001;Lenaerts et al., 2019). In a warming climate, however, this contribution could change sign to a net mass loss, accelerating the ongoing shrinking of ice sheets (Cullen et al., 2014). Such future changes, influencing SMB and energy balance, can only be predicted by climate models and rely on an accurate parameterization of humidity fluxes. Recent studies indicate, however, that climate models seem to systematically underestimate turbulent surface fluxes over ice surfaces (Box & Steffen, 2001;Fettweis, 2007;Noël et al., 2015).
Stable isotopologues of water act as natural tracers within the water cycle and provide an independent constraint on the parameterization of the exchange between the snow surface and the atmosphere. Each time a water parcel undergoes a phase change process, its characteristic isotopic fingerprint is modified. The isotopic composition of a water parcel can therefore be interpreted as an integrated signal of its path in the water cycle (Galewsky et al., 2016). As a consequence, the isotopic composition of, for example, the atmospheric water vapor can be used as an observable quantity to evaluate isotopologue-enabled model performance (e.g., Gryazin et al., 2014;Steen-Larsen et al., 2017). In isotope-enabled models, the individual water isotopologues are simulated independently. This way, difficult to observe processes like atmospheric entrainment and other mixing processes can be constrained through observations of water isotopologues (e.g., Benetti et al., 2018;Benson et al., 1994). Using multiple isotopologue-enabled General Circulation Models, Steen-Larsen et al. (2017) documented biases in the simulated water vapor isotopic composition above the GrIS despite only small differences between the observed and simulated specific humidity. They argued that this could indicate that the models did not accurately capture the moisture sources for the vapor above the ice sheet and hence illustrated how water isotopologue observations provide additional information on the modeling of the atmospheric hydrological cycle.
To identify the relevant processes responsible for potential model biases using water isotopic composition, the logical next step is to investigate the fluxes connecting the reservoirs with respect to their isotopic fingerprint. Such model verification requires accurate in-situ observations of the water isotopologue fluxes between the surface and the atmosphere. Those can be obtained either indirectly, by the Flux-Gradient (FG) method (Yakir & Wang, 1996) or the Keeling-Plot (KP) method (Keeling, 1958), or directly, by eddy covariance (EC) measurements (Swinbank, 1951). Based on the typical time resolution of the measurement systems for stable water isotopologues in the order of seconds, the FG method would be the natural choice. Measurements of CO 2 and H 2 O isotopologue fluxes using the FG method have been used to quantify ecosystem activity in areas with moderate humidity levels (Griffis et al., 2005;Lee et al., 2007Lee et al., , 2009Welp et al., 2008). Good et al. (2014) and Wei et al. (2015) used the isotopic signal of the humidity flux to differentiate between evapo-and transpiration using the KP method. Both approaches are based on the assumption that turbulent fluxes can be interpreted as gradient induced diffusion and are the foundation for flux parameterizations in models. The quality of these flux estimates are based on measurements of high precision and assumptions of source and sink homogeneity and the absence of counter-gradient transport (Griffis et al., 2008). However, those bulk methods experience inherent shortcomings, in particular when it comes to polar areas, characterized by low atmospheric humidity content and predominantly stable conditions (Anderson & Neff, 2008). The established gradient-flux relationships are based on the Monin-Obukhov similarity theory (MOST), but the underlying assumptions for the use of this approach may be frequently violated under stable conditions (Grachev et al., 2013;Mahrt, 2008;Sorbjan & Grachev, 2010). This issue can be expected to become more severe with increasing atmospheric stability in the surface layer. In addition, the gradient measurement strategy requires simultaneous measurements at different levels, and thus, optimally two analyzers that have to be carefully calibrated and maintained to provide the corresponding measurements with sufficient precision.
To overcome the inherent drawbacks of the gradient dependent approaches, the EC method is often applied. The EC technique is widely accepted as a direct and robust method for measuring latent and sensible heat fluxes (e.g., Foken et al., 2012;Lee et al., 2005) and has previously also been used for measurements of trace gas exchange (e.g., Baldocchi, 2014;Wohlfahrt et al., 2009) and isotopologue fluxes (e.g., Braden-Behrens et al., 2019;Griffis et al., 2008Griffis et al., , 2010. The EC method is based on the understanding of the eddy energy cascade which states that transport happens over a wide range of frequencies (Kolmogorov, 1991). This implies the demanding requirement of high precision and high temporal resolution measurements (Foken et al., 2012). Measurement systems, which do not meet the sampling requirements due to instrument WAHL ET AL.
limitations and set-up, have to rely on cospectral frequency loss corrections to obtain accurate flux measurements (Moore, 1986). Depending on the system set-up, both high and low-frequency corrections, might be required. Such corrections of measured covariances can be done by applying empirically derived models (Kaimal et al., 1972), utilizing analytical approaches (Massman, 2000;Spank & Bernhofer, 2008) or using on-site reference measurements for comparison, as outlined by Ferrara et al. (2012).
To establish isotopologue-enabled climate models as tools to improve our understanding of climate dynamics and climate projections, they need to have accurate process parameterizations implemented. Expanding on the limited numbers of stable water isotopologue flux measurements and adding observations, especially in the polar areas, is therefore a necessity to contribute to isotopologue-enabled climate model advancement. Unfortunately, this adds the complications emanating from low humidity and low temperature polar environments to the inherent challenges of scalar flux measurements. Hence, there is a need for a robust method to obtain stable water isotopologue fluxes.
This study presents a method for measuring stable water isotopologue fluxes in polar areas by combining high-frequency EC measurements with ∼1 Hz water vapor isotopologue cavity ring down spectroscopy (CRDS) measurements. CRDS instruments have previously demonstrated a good performance of continuous, high precision in-situ water vapor isotopologue measurements in high latitude locations (Berkelhammer et al., 2016;Bonne et al., 2014;Casado et al., 2016;Steen-Larsen et al., 2013). To assure robust isotopologue flux estimates in such environments, a stringent quality control based on integrated one-sided cospectra (ogive) analysis is outlined. Measurements obtained with this method will help to improve the process of understanding throughout the water cycle, whenever isotopologues are exchanged. We demonstrate the use of this method by presenting stable water isotopologue flux measurements between the snow surface and the atmosphere at the East Greenland Ice-Core Project (EastGRIP) site in the summer of 2019. Our results show that the isotopic composition of the humidity flux is dependent on vapor and snow isotopic composition, and that fractionation processes take place during sublimation.

Measurement Site
The measurement site was located at 75° 37′ 47′′ N, 35° 59′ 22′′ W in North-East Greenland at 2,700 m altitude on the GrIS. Measurements were conducted as part of the surface science division of the ongoing international deep ice core drilling project EastGRIP. Turbulent fluxes of heat, moisture, and water isotopologues were sampled continuously throughout the summer of 2019 from 25 May to 29 July. The local time (LT) corresponds to the UTC-2 timezone. Instruments for this study were installed on a tower in a dedicated clean snow area, with other camp and housing facilities downwind, to ensure undisturbed measurements ( Figure 1).
Radiation and relative humidity data were obtained from a continuously operating automatic weather station (AWS) from the Program for the Monitoring of the Greenland Ice Sheet (PROMICE) (Ahlstrøm et al., 2008;Fausto & van As, 2019), installed roughly 1 km SE from the flux measurement site. The synoptic conditions throughout the measurement campaign were characterized by a prevailing and rather stable wind direction of 285°N (±1  = 53°), a mean 2 m air temperature of −9.9 °C (max = 1.3 °C, min = −25.4 °C), low average specific humidity of 1.9 g/kg (max = 4.7 g/kg, min = 0.4 g/kg) and a distinct diurnal cycle visible in all relevant meteorological parameters on clear sky days.

Instruments
The instrumentation on the measuring tower used for monitoring isotopologue fluxes directly, consisted of a combination of two instruments. An IRGASON 2-in-1 sonic anemometer and open-path mid-InfraRed absorption Gas Analyzer (IRGA) by Campbell Scientific measured high-frequency three-dimensional wind and humidity and an optical Cavity Ring-Down laser absorption Spectrometer (CRDS), type L2140-i from Picarro measured stable water isotopologues in the vapor. The relative concentrations of naturally occurring stable water isotopologues 16 DH O (also D), 10.1029/2020JD034400 3 of 24 together with specific humidity by the CRDS. Specific humidity was reported in ppm v. Since the rare and heavy water isotopologues occur in small amounts compared to the abundant species ( 16 2 H O), the rare species were reported in ‰ following the standard  -notation after H. Craig (1961): R is the ratio of concentrations between the rare(*) and abundant isotopologue and * VSMOW R is the ratio of the reference Vienna Standard Mean Ocean Water (VSMOW) (de Laeter et al., 2003).
The IRGA was installed at an initial measurement height of 2.15 m above the snow (2 m-level), facing the prevailing wind direction. There was no significant change in the snow height observed throughout the season and all calculations were made using the initial height throughout the campaign. The IRGA measured 3D-wind components   ( , , ) U u v w and absolute humidity  w in the same measurement sample-volume. The CRDS inlet for vapor probing was mounted at 2.15 m, receiving the air sample just as it left the IRGA sampling volume with an offset of 10 cm in the prevailing wind direction (Figure 2). Three additional vapor inlets were installed on the tower at 0.47, 0.96, and 7.13 m height, referred to as the 0.5 m-, 1 m-and 7 m-level. The four inlets merged into one tube at the base of the tower using three three-way solenoid valves (EVT317-6DZ-02F-Q). From this junction, the air sample was pumped through a 15 m heated copper tube (OD: 1/4") using a KNF pump (N811 KN.18) with a flow rate of 10 L/min to the CRDS, which was kept in a tent at a constant temperature of  15 C. The overall measurement set-up is schematically presented in Figure 2.

Sampling Strategy
The IRGA measured wind, humidity, and sonic temperature continuously at a sampling frequency ( s f ) of 20 Hz. The raw data were stored on a CR6 Campbell Scientific Datalogger and the data were monitored using the Campbell Scientific program LoggerNet.
WAHL ET AL.   Vapor was pumped from four levels using the Picarro's valve sequencer software in a repeated one hour cycle of 10 min at levels 0.5, 1, and 7 m, and 30 min at the 2 m-level together with the IRGA high-frequency measurements. This resulted in a total of 1,350 30 min-periods for eddy covariance flux computation at the 2 m-level, hereafter called EC-periods. The CRDS ran in the 17 O-sampling-mode and vapor data of humidity and isotopologue concentration were generated at a s f of ∼1 Hz.
For comparison with the snow surface isotopic composition, half centimeter surface snow samples were sampled daily from 10 locations along a 100 m long transect. The samples were shipped frozen to the stable water isotope laboratory at Alfred-Wegener-Institut, Bremerhaven, Germany and analyzed on a CRDS as one consolidated sample following the Van Geldern and Barth (2012) measurement protocol.

Calibration and Uncertainty
The in-situ CRDS measurements of the isotopic composition of the vapor were calibrated against the VSMOW-SLAP scale using secondary laboratory standards fully spanning the range of observed water vapor isotopologue values provided by the Institute of Arctic and Alpine Research, University of Colorado (see Table S1). The CRDS instrument was calibrated once every day by providing a stream of vapor from individual standards using a custom made calibration system. As the CRDS instrument was observed to be very stable, the measurements were combined to one set of VSMOW-SLAP calibration curves for  18 O and  D to reduce the influence of the measurement noise from the day-to-day standard measurements. To account for the expected humidity dependency on the isotopologue measurements, humidity-isotopologue response curve calibration runs were performed several times during the campaign to test for stability . The humidity-isotopologue response curve was established by creating a stream of water vapor with a constant isotopic composition at different humidity levels . No drift in the humidity-isotopologue response curve was observed during the campaign, but the humidity-isotopologue response curve was dependent on the isotopic composition of the individual standards as previously observed by Weng et al. (2020) (Figures S6-S9). Linear interpolation in the isotopologue space between the humidity-isotopologue response curves was applied to correct for this dependency. The uncertainties on the CRDS vapor isotopologue measurements were assumed to be 0. To calibrate all available instruments for humidity, the AWS humidity data quality was tested and confirmed to be accurate by comparing the measurements with observations of an additional on-site CRDS analyzer. This analyzer had been calibrated against a custom made dew point generator spanning saturation vapor pressure levels from −40°C to 5°C. The average of the bias between the humidity measurements of these two was less than 0.001 g/kg for the specific humidity range of 0.4-1.6 g/kg. The bias in higher humidity conditions is assumed to be of the same order of magnitude. The AWS humidity measurements were then used as reference for CRDS and IRGA humidity post-season calibration.
Flux error computation required the specification of uncertainties of the individual parameters. The precision statements of the manufacturer for the IRGA were 0.006 mmol/mol on the humidity readings and 0.5 mm/s on the vertical velocity measurements. We defined the uncertainty of averaged vertical wind data to be the standard deviation of the mean over the averaging period.
The manufacturer does not give a precision value for the CRDS humidity readings, so the CRDS precision was estimated using controlled humidity calibration pulses. The obtained precision values for the isotopologues, as well as the humidity measurements, are conservative estimates of the actual instrument uncertainty, since the vapor stream, provided by the calibration system, was variable in itself. Furthermore, it is stated by the manufacturer, that the CRDS precision is humidity dependent. In agreement with this, we inferred a 0.2% uncertainty on the isotopologue as well as the humidity readings of the CRDS.

Flux Calculation Method
In the following section, we describe the protocol for calculating kinematic isotopologue fluxes ( kin

Preconditioning of Data
All 1350 EC-periods, in which CRDS and IRGA measured at the same height, were evaluated for isotopologue flux calculations. Data of the first and last 2.5 min were excluded from the beginning and end of each EC-period to avoid influences from the valve switching and allow for synchronization between instruments. The CRDS data set was interpolated to homogeneous 1 Hz-frequency. Humidity and isotopologue concentrations were converted to mass mixing ratios ( 18 , , CRDS D O q q q ) with units g/kg. Absolute humidity measurements of the IRGA (  w) were likewise converted into mass mixing ratio ( IRGA q ). The 20 Hz IRGA data were despiked following the procedure of Schmid et al. (2000) and periods with more than 1% of missing data were excluded (319 EC-periods discarded). Individual 25 min wind data sets were rotated into the mean wind direction using the double-rotation method, setting the mean w and v wind vectors to zero ). The despiked and rotated IRGA time series of w and IRGA q were downsampled to averaged 1 s-resolution for isotopologue flux computation.

Humidity Flux Computation
The IRGA humidity flux ( kin q IRGA F ) was calculated from the IRGA's mixing ratio and vertical wind data. Covariances (w q IRGA   ) were computed for a 25 min integration time ( ) after linear detrending. We evaluated shorter integration times, but found no significant improvement in the quality of the flux estimates. Integrated one-sided cospectra, called ogives (Og) (Desjardins et al., 1989), were computed from w and IRGA q for all despiked EC-periods following: 10.1029/2020JD034400 6 of 24 This presentation method of fluxes has two advantages. It illustrates how the total flux energy is distributed over the frequency range and it provides a method for evaluating the quality of the flux estimate with a given  . The endpoint of the ogive (i.e., the value at lowest frequency N f ) is equal to the covariance w q IRGA   with units in g m/(kg s) (Mann & Lenschow, 1994).
Frequency corrections applied to w q IRGA   comprised of low-frequency-loss corrections originating from linear detrending and high-frequency loss corrections due to path-averaging (Massman, 2000;Van Dijk, 2002;van Dijk et al., 2004). Webb-Pearman-Leuning corrections due to buoyancy effects (Webb et al., 1980) were not necessary, since fluxes were computed from the mixing ratio data directly (Businger, 1986). We used the established cospectral models of Kaimal et al. (1972)

Isotopologue Flux Computation
Due to the lower sampling resolution of the CRDS and likely additional tube attenuation effects, frequency corrections of covariances, computed from CRDS, were necessary. In this context, the kinematic humidity flux of the IRGA and its ogive is considered to represent the total humidity flux of each EC-period. It is used to compute the correction coefficient ( Og cc ) to correct calculated isotopologue fluxes. The IRGA humidity flux measurements themselves could be verified by good agreement with two other independent on-site EC-systems operated during the field campaign.
To account for expected lag times between the IRGA and CRDS data, 1s IRGA q and CRDS q were lag time-cross-correlated. The identified lag showed to be nonuniform in time throughout the season (Figure S1), which was due to a drift in time of the CRDS software system. Hence, respective lag times yielding maximum correlation between 1 s-IRGA q and CRDS q time series were used to shift the CRDS data accordingly. EC-periods with lag times exceeding the values within the observed trend were excluded (7 EC-periods discarded).
18 and corresponding ogives were calculated using the 1 s vertical wind data from the IRGA. The ratio between the end values of the humidity IRGA and CRDS ogives yields the Og cc as: All covariances computed from the CRDS data were corrected with the interval-specific Og cc . Finally, the frequency correction factor (FCorr) established for the humidity flux, was likewise applied to the CRDS isotopologue fluxes.
Knowing both, the total humidity and the isotopologue flux, one can calculate the isotopic composition of the humidity flux ( F ) as, which can be converted to  -notation as introduced in Equation 1.

Alternative Methods to Calculate Isotopic Composition of the Humidity Flux
In order to compare our estimates of the isotopic composition of the humidity flux with estimates from methods like the Keeling-Plot or Flux-Gradient approach, we recalculate  F with information from all four levels of the measurement tower. It is important to notice that neither of those two methods as presented below is suitable to estimate total isotopologue surface fluxes.

Flux-Gradient Method
For the FG approach, we made use of the 0.5 m-, 1 m-and 7 m-level 10 min measurements (Section 2.3) together with the 2 m measurements that were used in the EC calculations. Averages were computed for the three 10 min levels cutting the first 3 min after each valve switch, resulting in one measurement point at each level per hour. Likewise, 3 min were cut from the beginning of the 2 m measurement period. The remaining data were split into three subsets and averaged, yielding three sample points per hour at the 2 m-level. Subsequently, measurements at each level were linearly interpolated to the mid-point of each EC-period. Assuming equal turbulent diffusivities for all three isotopologues, the isotopic composition of the humidity flux can then be calculated from the ratio of heavy isotopologue gradient to light isotopologue gradient. Making use of information from all four levels, we calculated * ( ) F R FG as slope of the linear regression of CRDS q at each level against * q at each level using the following relationship (Good et al., 2012): WAHL ET AL.

Keeling-Plot Method
For the KP estimates of  F , we assumed a steady-state within each EC-period and the half hour that followed. We averaged the last 5 min of each level yielding four measurement points in 35 min. The KP approach is essentially a two end-member mixing model with one mixing source being the free atmosphere and the other mixing source being the humidity flux coming from the ground. Measurements at several levels in the atmospheric boundary layer then fall on the mixing line of those end-members. Consequently, and without knowing the actual value of the free atmosphere, we can estimate  * ( ) F KP by linear regression. Therefore, we establish the linear relationship between the reciprocal of CRDS q and  * and extrapolate to  The standard error of the regression intercept was taken as uncertainty estimate for the    KP F value.

Quality Assurance-Correction Coefficient Analysis
Calculation of fluxes using the EC method requires the check for compliance with the eddy covariance flux theory. The most important assumptions underlying the EC theory are homogeneity of the terrain, the absence of divergence or convergence and stationarity within the integration period (Foken & Wichura, 1996). Our measurement site satisfied the first two requirements, being a vast, open and homogeneous flat surface in the direction of the wind. Stationarity within the  was assured through the IRGA mixing ratio covariance (w q IRGA   ) dispersion of <30% following Foken and Wichura (1996), and periods with higher divergence were discarded (138 EC-periods discarded). Furthermore, we investigated possible measurement errors introduced by the measuring system by performing an extensive ogive analysis: First, CRDS-ogives of all three isotopologues were tested for their minimum responding frequency ( res f ) and where  norm denotes any normalized covariance from the CRDS. This was done to assure that Og cc can be used to correct fluxes of all three isotopologues.
Second, convergence analysis was performed on the ogives of the IRGA q -flux and the D q -flux as representative for the CRDS fluxes, following a similar approach as Foken et al. (2006). In contrast to continuous EC measurements at one height, our intermittent data set could not satisfy the originally proposed evaluation criteria, since a maximum analysis period of six times  would be required. Hence, a quality criterion based on the ogive shape was developed. For an adequate  , the ogive of a flux measurement will converge.
Consequently, if a log-linear fit through the low-frequency end of the ogive had a slope close to zero, that is, it converged, and this EC-period was accepted. The fit was calculated for ogive values with frequencies <4.5 * 10 −3 Hz using the least squares method. Ogives were characterized according to absolute values of the normalized slope ( norm m ) and normalized root mean square error ( norm RMSE ) of the fit as defined in Table 1. EC-periods with nonconvergent IRGA q -ogives were excluded from further flux calculation (9 EC-periods discarded). The   w D -ogives were classified as defined in Table 1. In QC2 EC-periods, the absolute value of the computed fluxes was slightly underestimated, as the ogive either had not reached its asymptote yet or had passed the extreme value already.

Isotopologue Flux Error Estimation
Conventionally, errors on EC flux measurements are either inferred from multiple instruments measuring the same flux or measuring during similar conditions, that is, comparing several days (Good et al., 2012). A more statistical approach was developed by Mann and Lenschow (1994) (ML) in which the integral length scale (  f ) and the correlation coefficient (r) between the time series of the vertical wind (w) and the scalar of interest (s) are considered.  f can be estimated by the ratio of measurement height to averaged wind speed within the  . The standard deviation of the flux estimation after ML then reads: with N being the number of samples in  and      / ( ) w s r w s . Finkelstein and Sims (2001) (FS) estimated the error on the flux obtained from EC measurements by calculating the variance of the covariance. Following this approach and adopting the notation from Fuller (1976), the standard deviation of the covariance between w and s is defined as: where the time lag (p) is defined by m = N/2 and is added to either w or s in the auto-(     / w w s s ) and cross-(     / w s s w ) covariances. For comparison, we calculated the relative error on the isotopologue flux values using both Equations 10 and 11. Since the calculated isotopologue fluxes are corrected using observed differences in humidity ogives between the two instruments ( Og cc ), the final flux isotopic value is essentially a product of three individually computed covariances. We therefore propagated the relative errors on the individual covariances, assuming they are not independent, to yield a final error estimation of the corrected isotopologue flux. However, this method produced rather high uncertainties that originate mainly from the uncertainty of the correction coefficient.
In order to differentiate errors stemming from instrument noise and sampling errors of the covariances (Businger, 1986), we set up an empirical error estimation using the Monte Carlo (MC) method. We carried out 1000 MC simulations for each EC-period that passed the quality control by adding random normally distributed noise based on the precision estimations of the respective instrument on the individual time series. The error on the final isotopologue flux value due to instrument noise was obtained from the standard deviation of the MC generated distributions. All three approaches were used to compute a relative error of the final isotopologue flux measurements and to separate instrument noise from stochastic (sampling) noise.

Measuring System Characteristics
We evaluated whether our sampling set-up resulted in unequal sampling of different water isotopologues due to potential tubing effects. For this, the minimum responding frequency ( res f , Equation 9), defining the noise and signal range, was calculated. The differences in minimum responding frequency between the total humidity and the heavy isotopologues were evenly centered around zero ( Figure S2). No systematic bias between the different water isotopologues was observed and the distribution around zero indicated randomness. Based on this analysis, we concluded that the isotopologues were not affected differently by the system. Hence, all measured isotopologue fluxes for a given interval could be corrected accurately with the calculated  of the individual water isotopologues measured with the CRDS. The averaged IRGA powerspectrum follows the expected slope for the inertial subrange. In contrast, the averaged CRDS isotopologue powerspectra show a dampening at high frequencies and a minimum responding frequency of ∼0.1 Hz, with higher frequencies being influenced by noise.

Ogive Correction Coefficient
For obtaining the Og cc , the humidity mixing ratio ogives of CRDS and IRGA of each individual EC-period were compared and evaluated according to Equation 5. The value of Og cc is directly linked to the shape of the ogives in the respective EC-period, and therefore the atmospheric stability. To demonstrate this dependency, ogives of example EC-periods are shown in Figure 5. The ogives of the CRDS humidity flux only start increasing at a frequency of ∼0.1 Hz, which corresponds to the signal-noise threshold as observed in Figure 4. In contrast, the ogives of the IRGA, being an open-path instrument and measuring at a higher frequency, start increasing at a frequency of ∼5 Hz.

Dependency on the Atmospheric Stability and Friction Velocity
Turbulence is created by the interplay of two influencing mechanisms: buoyancy and wind shear. The two corresponding key variables are the buoyancy flux (   w T ) and the friction velocity ( * u ). The friction velocity is always reinforcing a turbulent transport and increases with increasing wind speed. Depending on the WAHL ET AL.  sign of   w T , buoyancy can enforce or counteract turbulence. The ratio of shear and buoyant turbulent forcing is defined as the Monin-Obukhov length L as developed by Obukhov (1946). Taking into account the measurement height z = 2.15 m, the atmospheric stability can then be described by the stability parameter  = z/L. Unstable and stable atmospheric conditions are characterized by a negative or a positive stability parameter, respectively.
We used the accepted EC-periods (Table 1) to investigate the dependency of our method's correction coefficient ( Og cc ) on  and * u . We limited our analysis to a stability range of −0.3 <  < 0.06, which includes weakly stable regimes but excludes very stable conditions of  > 0.06 (Mahrt, 1998). We concentrate on this stability range, as transport in very stable conditions is likely not predominantly governed by the two turbulence parameters  and * u as further discussed in section 5.1.1. To analyze the dependency of Og cc on * u , we analyze the stability range of − 0.3 <  < 0.01, as the stability influence would otherwise mask the friction velocity dependency. The reader is referred to the supporting information to see the same analysis including very stable conditions and the dependency on wind speed ( Figures S3 and S4).
The Og cc was found to increase exponentially with increasing stability as shown in Figure 6a). When focusing on the near-neutral and unstable regimes (−0.3 <  < 0.01), a dependency of Og cc on the friction velocity becomes apparent. Values of Og cc increased linearly with the friction velocity (Figure 6b).

Quality of Data
During the analysis of the low-frequency shape of the   w D -ogives (Section 3.5) 578 of the evaluated EC-periods were accepted; 371 as QC1 and 207 as QC2 (Table 1). A dependency on stability was observed. In stable conditions, the measured   w D -ogives did not develop a pronounced S-shape which would indicate that both, the high and low frequency fluctuations were appropriately covered by the chosen averaging period and the instrument sampling frequency. Therefore, those ogives had small end values close to zero resulting in corresponding high

Flux Error Analysis
The three different error calculation methods resulted in three clearly separated distributions for the relative error of the final isotopologue flux values (Figure 7). All three error distributions are skewed to the right but their median is located at different values for the error sum of 5%, 29%, and 56% for the MC, ML, and FS approach, respectively. Generally, high relative error values of the ML and MC methods could predominantly be associated with a stably stratified atmosphere and low absolute flux values, whereas values from the FS method could not readily be assigned to one regime only. The errors obtained from the MC method were an order of magnitude smaller than the errors calculated by the FS method. This was to be expected, as the MC approach only takes into account instrument noise and omits stochastic error contributions, that is, the number of independent observations. Both the ML and FS method are estimates of the total uncertainty of the flux measurements.

Computed Isotopologue Flux Values
Three clear sky periods (one each in May, June, and July) were chosen to demonstrate and compare the behavior of isotopologue surface fluxes during ideal summer conditions in the interior of the GrIS (Figure 8). In all three periods, the averaged net latent heat flux (LE, Figure 8c) was slightly positive, that is, sublimation outweighed deposition and the humidity flux contributed negatively to the surface mass balance. Likewise, the net surface flux of isotopologues was directed away from the surface and the snowpack lost heavy isotopologues to the atmosphere through sublimation. In general, the isotopologue fluxes co-varied with LE. Between May and June, the amplitude of the 18 O-isotopologue surface flux increased from an averaged daily maximum value of 0.5 10 −5 g/kg to 1.5 10 −5 g/kg (Figure 8a). The measured D-isotopologue flux was of an order of magnitude smaller than the measured 18 O-isotopologue flux but showed the same pattern ( Figure 8b). This increase in mass exchange was also observed in the sublimation rate (calculated from LE assuming a latent heat of sublimation of 2.834 10 6 J/kg) with averaged net daily values during the June period (0.14 mm water equivalent [w.eq.]) being five times larger compared to the averaged 0.03 mm w.eq. that sublimated during the May period. During July, the net sublimation rate decreased again to 0.08 mm w.eq. The mass exchange increase between May and June is accompanied by higher temperatures (Figure 8e) and increasing humidity. Absolute humidity observations (not shown) recorded an increase in water vapor content throughout the season from values of around 1 g/m 3 in May to 3 g/m 3 in July. The air temperature and relative humidity were near inversely proportional and on diurnal time scales driven by the net radiation ( Figure 8h). Consequently, and characteristic for polar areas, the air had a general high relative humidity level (mean: 86%) despite the low absolute humidity. Saturation (with respect to ice) was repeatedly reached during night time hours in the May period. However, it was not reached in the June period and only registered around local midnight for the July period, resulting in longer periods of sublimation conditions (Figure 8f). A weak diurnal cycle signal in the wind speed was observed but no considerable change in wind conditions between the three periods was recorded (Figure 8g). In all three periods, the average net sensible heat flux (H, Figure 8d)   was supported by snow temperature measurements (not shown). Strong negative daytime, but weak positive nighttime, temperature gradients between snowpack and air temperature resulted in high H daytime values and only slightly negative daily mean H values of −0.7 W/m 2 in May, when the air temperatures were still low. Further into the summer, a decreased net H flux showed stronger nighttime snow-air temperature differences during the clear-sky periods and a corresponding amplification of the energy uptake of the snow with a daily mean H of −5.8 W/m 2 and −5.4 W/m 2 in June and July, respectively.

Computed Isotopic Composition of the Humidity Flux
Although, in general, the isotopologue fluxes co-varied with the humidity flux, the relative abundance of isotopologues in the humidity flux changed depending on the underlying process. The identification of influencing parameters required the knowledge of the isotopic composition of the humidity flux. Therefore,  F was calculated with the eddy-covariance method using the 1 Hz CRDS measurements following Equation 6. For comparison,  F was also calculated assuming the Flux-Gradient theory ( F (FG), Equation 7) and following the Keeling-Plot approach ( F (KP), Equation 8). A detailed quality control procedure for the EC method was outlined in section 3.5 and summarized in Table 1. For the KP and FG approach, the estimated uncertainty of the  F value was dependent on the 2 R value of the linear regression and on the humidity difference between the observation levels (Δq) within the calculation period, as previously observed by Santos et al. (2012). High uncertainties were observed in the KP calculations when 2 R < 0.9. We therefore used WAHL ET AL.  Hourly surface flux measurements of latent (c) and sensible heat (d), ambient 2 m air temperature (e) and 2 m wind speed (g) are reported from the flux tower, whereas 2.5 m relative humidity (f) and net surface radiation (h) were measured at the PROMICE station. The sign convention is positive for upward fluxes. Note the discontinuous time axis. 2 R > 0.9 as the quality criterion for the KP estimates. Concerning the second criterion, Δq of the FG data was always smaller than Δq of the KP data. This is because we interpolated the humidity levels to the same point in time for the FG method and the corresponding Δq therefore does not include a time component unlike Δq of the KP method. We inferred a quality criterion of Δq > 50 mg/kg −1 for the FG method from the uncertainty analysis. Of the 578 available QC1 & QC2 EC-periods, 36% fulfilled the two criteria for  F D and 28% for  18 O F ( Figure S6).
The result of the comparison between  F (EC),  F (FG), and  F (KP) for sublimation conditions (LE > 0) is shown in Figure 9. Only the reduced data set, that passed the quality criteria, is shown. In general, the FG estimates (Figures 9b and 9d) had a higher correlation (r > 0.70) with EC results than KP estimates ( r > 0.65) (Figures 9a and 9c). Furthermore,  F D (Figures 9c and 9d) showed a slightly better agreement between methods than  18 O F (Figures 9a and 9b).
During deposition conditions (LE < 0),  F between the three methods compared less well with values not falling on the one-to-one line (Figures S7-S8). The correlation coefficients were smaller than 0.55 in both comparisons and for both isotopologues.
In the course of the season, isotopic enrichment of the sublimation flux was recorded. During sublimation conditions within the three example periods, the averaged isotopic composition from the EC method was Shown are the share of KP and FG estimates that passed the quality control (28% for  18 F and 36% for  F D) and the corresponding correlation coefficients (r). The estimates were matched on hourly resolution and the 1:1 line is shown in gray.
0‰ for  18 O and from −500 to 0‰ for  D. To account for the expected noise in the signal, the data was bin averaged ( 18O bin size: 1‰,  D bin size: 5‰) and the error bars indicate the standard deviation in each bin. The correlation between snow and flux isotopic composition was higher for  18 O (r = 0.81) than for  D (r = 0.65). Weighted linear regressions were performed for both isotopologue species and were significant. The corresponding equations and p-values are given in Figure 10. The sublimation flux was isotopically depleted compared to the snow surface isotopologue content. For comparison, we calculated the expected isotopic composition of vapor under the assumption of ice-vapor equilibrium fractionation. We used the observed mean surface temperature during sublimation conditions (−9 °C) and fractionation factors of Majoube (1971), Merlivat and Nief (1967) and Ellehoj et al. (2013) for this calculation.
For deposition conditions, a significant dependency of  F on the isotopic composition of the 2 m vapor ( V ) was observed (Figures 10b and 10d). Contrary to sublimation conditions, the correlation coefficient was used the observed mean surface temperature during deposition conditions (−14°C) and equilibrium fractionation factors of Majoube (1971), Merlivat and Nief (1967) and Ellehoj et al. (2013) for the calculations.

System Performance
Our measurement system has proven to produce reliable flux measurements of all three water isotopologues (Section 4.1, Table 1). Despite small differences in molecular characteristics between the different isotopologues, we found no differentiation by the measurement system. This agrees with assumptions made by Good et al. (2012) following Griffis et al. (2010), who found similar behavior for

Correction Coefficients for the EC Method
The lower sampling resolution of the CRDS made corrections necessary to account for high-frequency flux contributions that could not be resolved by the instrument. Specifically, the CRDS′ low sampling resolution can be regarded as a low-pass-filter. As the CRDS′ sampling frequency is only a fraction of the IRGA's sampling frequency and tube attenuation likely dampened the high-frequency signal additionally, the calculated Og cc -values were high. However, this should be regarded as an inevitable characteristic, rather than a deficiency of the method. Despite high maximum values of Og cc up to 15 during moderately stable conditions, and even higher values for a strong atmospheric stability, the coefficients showed to be systematically dependent on turbulence parameters and can be explained as follows: Depending on the state of the atmospheric boundary layer, relatively more or less energy is transported at higher frequencies, corresponding to a higher or lower correction coefficient, respectively. The stronger the stable stratification, that is, the higher the  -value, the smaller the dominant eddies with corresponding higher frequencies. This means that the cospectral peak shifts to higher frequencies with increasing atmospheric stability. Due to the lower sampling frequency, the CRDS system cannot resolve this progressively bigger high-frequency share of the flux and the correction factor therefore becomes larger.
The lowest Og cc -values were observed in unstable stratification. Here, the system was able to resolve most of the lower frequencies, at which transport is occurring under these conditions. Although generally lower, the cospectral peak under unstable conditions likewise shifts to higher frequencies with increasing wind speeds. This explains the increasing Og cc with increasing friction velocities under unstable or near neutral conditions. In stable regimes, this dependency is masked by the pronounced effect of stratification on the cospectral peak (Massman, 2000), and the reason why stability regimes of  > 0.01 were excluded in the analysis. A linear dependency of correction factors on wind speed, was previously observed by Ferrara et al. (2012) for empirical and analytical frequency correction approaches of CH 4 -fluxes. However, it is important to note, that this behavior is intrinsic to closed-path analyzers like the CRDS. Massman (2000) indeed, found a slight decrease of correction factors with higher wind speeds for open-path instruments that suffer merely from high-pass-filters like block averaging.
Despite the reasonable justification for high Og cc -values during stable atmospheric conditions, it needs to be emphasized that stable stratification is likewise the biggest hindrance for this type of flux measurements. Most EC-periods that did not pass the ogive-quality test were periods during stable and very stable conditions. Previous studies have therefore often excluded stable stratification periods completely (Good et al., 2012). However, with our site of interest being predominantly stably stratified, we included stable conditions in our analysis and applied the same quality criteria to both, stable and unstable conditions. In doing so, we were able to identify reliable data to draw conclusions on the vapor isotopologue exchange through sublimation and deposition during polar day conditions, even in nonconvective regimes. Based on this, it needs to be mentioned that instruments with a higher temporal resolution will be necessary to measure isotopologue fluxes during polar night conditions. Strong atmospheric stability, as observed during polar night conditions, likely requires shorter integration times of seconds to minutes. Combined with the low resolution, closed-path analyzers, this limits the available frequency range, from which flux information can be gained. Furthermore, it needs to be kept in mind that the EC method is developed for direct WAHL ET AL.

10.1029/2020JD034400
17 of 24 measurements of turbulent transport. However, under a strong stable stratification, transport might not comply with the turbulence theory alone but might also include submeso-motion-like transport (Acevedo et al., 2016;Mahrt, 2009). The general discussion of transport and its measurement in stably stratified conditions lies, however, beyond the scope of this work.

Comparison of  F Between Methods
Our EC estimates of the isotopic composition of the flux compared well to the more established FG and KP methods. During sublimation conditions, when uncertainties of EC estimates are smallest, FG and EC showed a correlation of >0.7. Correlation coefficients between the EC and KP methods were slightly smaller which could be due to the timing offset between the 25 min EC-period and the subsequent 35 min KP period. The steady-state assumption, necessary to compare the estimates of the two methods within one hour, might not be valid which could explain the difference in the estimates. For a more direct comparison, simultaneous measurements at several heights with multiple instruments would be needed which was not feasible with the available system set-up. In order to bypass the time offset, we interpolated the measurements of all levels to the midpoint of the EC period for the FG approach. This way, we introduced the uncertainty of interpolation but satisfied the steady-state assumption. This could be the reason why the FG results agreed better with the EC results. We point out, that if the exact same data were used for FG and KP estimates, both methods would yield the same result as noted by Good et al. (2012).
The FG and KP pitfalls outlined earlier become apparent when looking at the number of available periods for the comparison. Over half of the periods that passed the EC ogive quality screening did not pass the KP R 2 and FG humidity range (q) quality screening. This alone might speak for the EC method. However, as we have only evaluated available EC-periods for the comparison, FG and KP methods might be eligible at times when the EC quality criteria were not met. Nevertheless, the EC is a direct method, which allows measurements of fluxes without relying on simplified conceptual theories. As the stable water isotopologue analyzer development proceeds, an increased sampling frequency combined with reduced memory effect and smoothing will allow for the use of the full potential of the EC method when measuring isotopologue fluxes. Besides the high sampling frequency, an additional instrument requirement is the very low instrument drift. One quality criterion for EC-period rejection was defined as a high noise level on the low-frequency end of the cospectrum (  0.1 norm RMSE ). Since the ogive end-value is the critical parameter for correction factor calculation, this value needs to be accurate. Mammarella et al. (2010) argue that the low-frequency noise stems from instrument drift. Hence, it is important that instruments used for this kind of flux measurements have minimum drift in their readings. As CRDS instruments have shown to have a high degree of stability despite the harsh polar conditions, we have thus used a CRDS analyzer in our study.

Error Estimation
We used three different methods to quantify and characterize the error on the EC isotopologue flux measurements. The MC method was used to inform about the random error contribution as it only takes into account instrument noise. When comparing the median of the MC distribution with the other two median values, it becomes obvious that instrument noise accounts for only a fraction of the total uncertainty. The MC median is several times smaller than the ML or FS median. Both ML and FS error estimations are methods, that calculate the variance of the covariance and therefore account for stochastic sampling errors (Businger, 1986).
Stochastic sampling errors are directly related to the amount of available independent observations. Averaged covariances need significantly more independent observations (i.e., longer integration times) to yield a specific low uncertainty compared to regular scalar averages. According to Businger (1986), this is due to the intermittent character of second-order moments which seem especially related to large eddies. The CRDS instrument predominantly measures these large eddies with low frequencies. It is therefore not surprising that we find such high error estimations for the covariances computed from CRDS time series. By computing the cross-and autocovariances using different lags, the FS method specifically quantifies the influence of the number of independent observations. The ML approach approximates, instead, the number of independent observations by the integral length scale. This is likewise a possible explanation why the FS method yields higher error values than the ML method. That said, total uncertainties obtained with both WAHL ET AL.
If we interpret the random noise contribution as a precision estimate of our method, it becomes clear why prominent diurnal cycle structures can be observed in our measurements despite the high uncertainty estimates of FS (∼50%) and ML (∼30%). The bigger share of the overall uncertainty would then be attributed to sampling errors, and the FS and ML uncertainties can be interpreted as the combined uncertainty of the flux values due to precision and accuracy. Since the system set-up was not modified during the campaign and all measurements were post-processed following the same protocol, we argue that this systematic sampling error does not change during the observation period and we can therefore compare measurement values from different periods within the observation period. The direct comparison of data from a second measurement system, however, has to be done with care, accepting potential differences within the range given by ML or FS uncertainty estimates.

The Surface Isotopologue Flux at EastGRIP
We concentrate here on three examples of clear sky conditions in the interior of the GrIS. Other clear sky period observations were similar to the ones shown. We chose these time periods for our focus as they represent idealized conditions when the diurnal cycle is the dominant mode of variability for the meteorological parameters and influences of more complex drivers, such as clouds and snowfall, are reduced. At the same time, clear sky conditions allow for an unstable regime of the lower atmosphere with a negative temperature gradient close to the surface. Such convective boundary layers have been observed previously on ice sheets (Anderson & Neff, 2008). Our method has shown to be more reliable during such unstable atmospheric stratification.
Those ideal clear sky conditions show a signal of isotopologue surface fluxes that follow the diurnal cycle of the net surface radiation but are limited in amplitude by the saturation water vapor deficit. With warming conditions, the saturation vapor pressure becomes higher and the air can hold more moisture, which likewise allows for higher latent heat and isotopologue flux values. Indeed, the isotopologue fluxes mirror the latent heat flux behavior and all fluxes increase after the May period. However, wind conditions do not change substantially which points toward a stronger buoyancy effect that enhances turbulence in the June and July period and can explain the increased fluxes.
In all three example periods, the net isotopologue fluxes and latent heat flux are directed away from the surface, which implies a net sublimation flux in accordance with observations from other locations in the accumulation zone of the GrIS (Box & Steffen, 2001;Cullen et al., 2014;van den Broeke et al., 2011), but contrary to some modeling studies (Box et al., 2006;Franco et al., 2013).
For identifying model biases, the isotopic composition of the water vapor flux is the parameter of interest. We find that the isotopic composition of the humidity flux during sublimation and deposition shows a dependency on the snow surface and vapor isotopic composition, respectively. The observed isotopic enrichment of the sublimation flux throughout the season can thus be explained by enrichment in the surface snow isotopic composition as new snow is deposited. Furthermore, recent process model studies (Madsen et al., 2019;Ritter et al., 2016) as well as laboratory experimental studies (Hughes et al., 2021) have suggested changes in snow isotopic composition due to isotopically fractionating exchange-processes. Our measurements of the isotopic composition of the sublimation flux support the hypothesis of sublimation induced fractionation. The calculated linear regression models fitted to observations do not follow a one-to-one line with the snow isotopic composition, which would be expected if sublimation did not cause fractionation as stated by Dansgaard et al. (1973). Instead, we observe a depleted isotopic composition in the humidity flux compared to the snow surface isotopic composition. The regression models for  D and  18 O are close to the theoretical equilibrium vapor lines that we calculated assuming established vapor-ice fractionation factors (Ellehoj et al., 2013;Majoube, 1971;Merlivat & Nief, 1967). This suggests that the isotopic composition of the sublimation flux can, to a first order, be calculated from the snow isotopic composition assuming equilibrium fractionation. Depending on the ratio between the sublimation and snowfall amount, the sublimation flux modifying the snow surface layer, might influence the snow isotopic signal that is later buried and stored in the snow column. Consequently, and determined by the meteorological conditions and pre-WAHL ET AL.

10.1029/2020JD034400
19 of 24 cipitation seasonality at ice core drilling sites, the fractionation surface processes could have implications for paleoclimate proxy interpretation of ice core records. It has to be investigated to which extend the sublimation signal overprints the precipitation climate signal and which meteorological parameters are critical for controlling the magnitude of shifts in the snow isotopic composition. Establishing an understanding of the exact mechanisms and influencing parameters will be the subject of future studies.
For deposition conditions, the isotopic composition of the humidity flux becomes relatively more depleted with the decreasing isotopic composition of the vapor. Under the assumption of steady-state, we expect that the isotopic composition of the deposition flux is essentially the same as the isotopic composition of the frost forming during deposition. Yet the regression models for both isotopologue species are not consistent with the expected equilibrium fractionation vapor-solid values. First, this is because we calculate the solid isotopic composition with regard to the 2 m-level although we expect a substantial gradient between the 2 m-level and the vapor very close to the surface. Casado et al. (2016) found steep gradients of the isotopic composition in the centimeter range near the surface suggesting a more depleted vapor close to the ground compared to the 2 m observations. This would in general lead to the observed  F being more depleted than what would be expected from equilibrium fractionation of the 2 m water vapor level. Second, differences between the slopes of the regression and equilibrium lines can be explained by influences of kinetic fractionation processes during deposition conditions. During these periods, as shown in Figure 8, the air is frequently saturated and, with respect to the surface temperature, at times supersaturated. Supersaturation, as shown by Jouzel and Merlivat (1984), introduces a humidity gradient at the condensation interface and therefore causes kinetic fractionation influences. Since supersaturation and the vapor isotopic composition very close to the ground are difficult to measure, parameterizations of  F are challenging to establish. The observed dependency between the 2 m vapor isotopic composition and  F could serve as an empirical benchmark in future surface exchange studies and for quantifying the role of supersaturation on kinetic fractionation.
From our observations, we can conclude that isotopic fractionation during sublimation and deposition alters the snow surface isotopic composition. Whether this atmosphere-snow exchange affects the long-term mean isotopic composition of the surface snow, and its implications for paleoclimate ice core proxy interpretations requires further analysis.

Conclusions
The combination of a high-frequency sonic anemometer with a high-precision slow sensor for isotopologue measurements allowed direct measurements of in-situ surface fluxes of stable water isotopologues using the eddy-covariance method. This has previously not been conducted in such low humidity conditions as found in the interior of the Greenland Ice Sheet. We make use of the available high-frequency cospectrum of a co-located humidity sensor to correct for the high-frequency-flux loss, which does not require any kind of similarity assumption. Our method allows for the quantification of the stable water isotopologue exchange between the snowpack and the atmosphere with the commonly known limitations of eddy-covariance flux estimations. In particular, strong atmospheric stability was identified as a hindrance to isotopologue flux measurements since correction factors of over 15 are needed in very stable conditions. In addition, calculated flux errors were found to be higher in stable conditions. This implies that the application of the presented method and instrumental set-up is most applicable during dominantly unstable to slightly stable conditions. When evaluating other possible instruments as a substitute to the slow sensor used in this study, we highlight the need for very low drift of the analyzer to minimize the low-frequency noise.
Surface isotopologue fluxes during clear sky conditions were readily observable and allowed a quantification of the isotopologue exchange at the surface. During polar day, the net flux of isotopologues was directed away from the surface, that is, sublimation outweighed deposition and the snowpack lost heavy isotopologues to the vapor. The isotopic composition of the sublimation flux showed to be dependent on the surface snow isotopic composition and was, when compared to the snow, isotopically depleted. This depletion can, to a first order, be approximated assuming equilibrium fractionation during sublimation. Furthermore, the observed isotopic composition of the humidity flux during deposition conditions was dependent on the vapor isotopic composition. Our findings imply that, depending on the strength of vapor-snow exchange, the surface fluxes can modify the precipitation-weighted climate signal in snow with implications for ice core WAHL ET AL.
10.1029/2020JD034400 20 of 24 records. Thus, the presented observations of the exchange of stable water isotopologues between the snow surface and lower atmosphere can serve as the basis for improved parameterizations of surface exchange processes in climate models and as an important component for snow isotopologue studies.