The influence of snow sublimation and meltwater evaporation on δ D of water vapor in the atmospheric boundary layer of central Europe

Post-depositional fractionation of stable water isotopes due to fractionating surface evaporation introduces uncertainty to various isotope applications such as the reconstruction of paleotemperatures, paleoaltimetry, and the investigation of groundwater formation. In this study, we investigate isotope fractionation at snow-covered moisture sources by combining 17 months of observations of isotope concentration ratios [HD16O] / [H16 2 O] in low-level water vapor in central Europe with a new Lagrangian isotope model. The isotope model is capable of reproducing variations of the observed isotope ratios with a correlation coefficient R of 0.82. Observations from 38 days were associated with cold snaps and moisture uptake in snow-covered regions. Deviations between modeled and measured isotope ratios during the cold snaps were related to differences in skin temperatures (Tskin). Analysis of Tskin provided by the Global Data Assimilation System (GDAS) of the NCEP implies the existence of two regimes of Tskin with different types of isotope fractionation during evaporation: a cold regime with Tskin < Tsubl,max = −7.7 C, which is dominated by non-fractionating sublimation of snow, and a warmer regime with Tsubl,max < Tskin < 0 C, which is dominated by fractionating evaporation of meltwater. Based on a sensitivity study, we assess an uncertainty range of the determined Tsubl,max of−11.9 to−2.9 C. The existence of the two fractionation regimes has important implications for the interpretation of isotope records from snow-covered regions as well as for a more realistic modeling of isotope fractionation at snow-covered moisture sources. For these reasons, more detailed experimental studies at snow-covered sites are needed to better constrain the Tsubl,max and to further investigate isotope fractionation in the two regimes.


Introduction
The hydrological cycle of the atmosphere is usually investigated by water vapor concentration measurements.A new dimension is opened by the analysis of H 2 O stable isotope ratios, which are modified by H 2 O phase changes during evaporation, cloud formation, and in-cloud physics.The main reason for this fractionation is the difference between vapor pressures of stable isotopologues such as H 16  2 O and HD 16 O, which results in a preferential condensation of HD 16 O.In consequence, not only specific humidity and dew point temperature but also the isotope concentration ratio R D = [HD 16 O]/ [H 16  2 O] -commonly referred to as δD = R D /R D,VSMOW − 1 with R D,VSMOW = 0.00031152 -decrease in a cooling and raining air mass.The resulting relation between condensation temperature and isotope ratios of water vapor or precipitation is the basis for a variety of applications.Water isotopes in ice cores are used for the high-resolution reconstruction of paleotemperatures (Dansgaard, 1964;Masson-Delmotte et al., 2005).The temperature-induced gradient of isotope ratios with altitude makes water isotopes in precipitation a proxy for paleotopography (Poage, 2001;Rowley et al., 2001;Blisniuk, 2005;Rowley and Garzione, 2007).Hydrological studies exploit the altitude dependence and a seasonality of isotope ratios in Published by Copernicus Publications on behalf of the European Geosciences Union.
precipitation to investigate groundwater formation (de Vries and Simmers, 2002).
Isotope fractionation during evaporation at the ground level post-depositionally modifies the isotope ratio of water at the surface and in the upper soil layers (Barnes and Allison, 1983).Because the various isotope applications rely on a close relation between the isotope ratio of precipitation and the isotope ratio of water at the surface, postdepositional effects increase the uncertainty of the isotope applications.Most problematic in that context are systematic post-depositional changes of isotope ratios into one direction, which most likely affect water reservoirs with long exposure to the atmosphere such as water intercepted in the snowpack.
Several studies report post-depositional enrichment of heavy isotopes in the snowpack.Such enrichment was observed over a wide range of temperatures and even at skin temperatures far below the freezing point.Epstein et al. (1965) observed enhanced δD values in depth-hoar layers in South Pole firn.To explain the observations, the authors suggested non-fractionating sublimation with part of the sublimated vapor fractionating recondensing and part of the sublimated vapor escaping the snow layer.Moser and Stichler (1974) analyzed surface layer snow from Switzerland (2450 m a.s.l.) during an 8-day fair weather period with air temperatures between about −5 and 0 • C.They observed a continuous increase of δD values in the aging surface layer snow, which they attributed to fractionating evaporation or sublimation of snow.Stichler et al. (2001) report a similar experiment from the Chilean Andes (5536 m a.s.l.) with air temperatures between about −12 and −5 • C.They observed an increase of δD values in the aging surface layer snow as well.They attributed this increase to kinetic fractionation during sublimation, caused by smaller coefficients of diffusion of the heavier isotopes.Based on a seasonal increase of δD values in a snowpack in northern Norway (∼ 900 m a.s.l.), Gurney and Lawrence (2004) suggested fractionating evaporation of meltwater and subsequent recrystallization of residual meltwater.Consistent with this, Lechler and Niemi (2011) report an altitude-dependent seasonal modification of δD values in an alpine snowpack, which they attributed to fractionating evaporation of meltwater during ablation.
A clear assignment of such observations to a specific type of interaction between snowpack and atmosphere or to certain meteorological conditions is difficult, as continuous time series of high-resolution isotope profiles of the snowpack are hard to obtain.Furthermore, processes such as meltwater percolation, diffusion, sublimation, and deposition within the snowpack additionally modify the profiles of isotope ratios, making an interpretation even more challenging.The current understanding of fractionation during the evaporation of snow is therefore highly uncertain and stretches from non-fractionating layer-by-layer sublimation of snow without any modification of isotope ratios in the snowpack (Ambach et al., 1968;Dansgaard, 1973;Friedman et al., 1991) to systematic enrichment of heavy isotopes in the snow in consequence of fractionating evaporation of meltwater (Gurney and Lawrence, 2004;Lechler and Niemi, 2011).This uncertainty complicates interpretation of isotope records from snow-covered regions, as the contribution of the various potential post-depositional effects may be different during climatologically different time periods.In addition, this uncertainty limits the isotope modeling of atmospheric moisture sources in snow-covered regions, as isotope-enabled models in general only consider one of the types of isotope fractionation for the sublimation of snow (e.g., Yoshimura et al., 2006;Werner et al., 2016).
In this context, continuous observations of isotope ratios of low-level water vapor may provide new insights by offering the opportunity to investigate fractionation during evaporation from the snowpack from a complementary point of view.A case study by Noone et al. (2013) investigates local variations of δD values in water vapor subsequent to a winter storm.Based on observations on a research tower, the study partitioned the different surface fluxes by assuming local evapotranspiration to consist of non-fractionating sublimation and fractionating evaporation of meltwater.
A promising way to extend the investigation of isotope fractionation during evaporation of snow to the remote moisture source regions is the combination of isotope observations with Lagrangian isotope modeling.A number of earlier studies applied Lagrangian isotope modeling along idealized climatological trajectories (Jouzel and Merlivat, 1984;Johnsen et al., 1989) or along individual back trajectories (Schlosser et al., 2004;Helsen et al., 2004Helsen et al., , 2005Helsen et al., , 2007;;Sodemann et al., 2008a;Pfahl and Wernli, 2009).As these studies mainly focused on polar or marine regions, the applied Lagrangian isotope models are not optimized to simulate continental evapotranspiration in the middle latitudes or evaporation of snow or meltwater at temperatures close to the freezing point.
In this paper, we present a time series with 17 months of continuous measurements of δD in low-level water vapor in central Europe.The measurements cover two winters, which were marked by a number of cold snaps and related snowfall.By combining the δD observations during the cold snaps with a new Lagrangian isotope model, we investigate isotope fractionation at snow-covered moisture sources.
In the following, we present the Lagrangian isotope model (Sect.2) and characterize the uncertainty of our δD measurements (Sect.3).In Sect.4, we relate variations of the δD to air temperature and specific moisture source regions during cold snaps.In Sect.5, δD values during cold snaps are analyzed with respect to isotope fractionation during surface evaporation in snow-covered regions.

Back trajectories
Isotope ratios to be analyzed in this paper were measured at a site near Karlsruhe in central Europe (49.10 • N, 8.44 • E; 110.4 m a.s.l.).Kinematic 5-day back trajectories from the site were calculated with the Hybrid Single Particle Lagrangian Integrated Trajectory model 4.0 (Draxler and Hess, 1998) with a time resolution of 1 h.Three-dimensional wind fields for trajectory calculation were derived from the Global Data Assimilation System (GDAS) of the National Centers for Environmental Prediction (NCEP) (Kanamitsu, 1989;Derber et al., 1991;Parrish and Derber, 1992).The GDAS data (V1.5) are available for every 3 h at a 1 × 1 • horizontal resolution and on 23 sigma pressure levels between 1000 and 20 hPa.To account for uncertainty of the back trajectories, we used trajectory ensembles.Each ensemble consists of nine trajectories, starting 30 m above ground level at the measurement site as well as 50 km N, NE, E, SE, S, SW, W, and NW from the site.Back trajectories were calculated for every 3 h of measurement time and were initialized at 00:00, 03:00, 06:00, 09:00, 12:00, 15:00, 18:00, and 21:00 UTC.

Lagrangian diagnostic of moisture sources
In order to identify major source regions of low-level water vapor in Karlsruhe, we applied a Lagrangian source region analysis similar to the method described by Sodemann et al. (2008b).The method traces air parcels along kinematic 5-day back trajectories and analyzes changes of specific humidity (q) in time intervals of 1 h.Such changes are possible due to the formation of precipitation (P ), evaporation from the ground (E), evaporation of falling rain, air mass mixing due to convection or small-scale turbulence, diffusion, and numerical errors.In addition, using GDAS wind fields with 3hourly time resolution for the calculation of back trajectories -which is a coarse time resolution compared to the GDAS time steps on the order of minutes -may cause deviations between the HYSPLIT back trajectories and the exact trajectories that air parcels followed in the GDAS.This, in turn, may result in artificial changes of q along the HYSPLIT back trajectories.For instance, HYSPLIT trajectories not fully capturing diurnal vertical movement in consequence of thermal expansion of the atmospheric boundary layer (ABL) in the GDAS may show artificial diurnal changes of q.Like Stohl and James (2004) and Sodemann et al. (2008b), we assume P and E to be the dominant processes and ignore the other effects.To avoid an overestimation of the formation of precipitation and moisture uptake in consequence of potentially artificial diurnal variations of q along the HYSPLIT trajectories, we smoothed q along the back trajectories with a 24 h rectangle function.Using these simplifications, the change of specific humidity ( q) per time step ( t) is q/ t = (E − P )/ t. (1) Assuming further that either E or P are dominating (James et al., 2004), the net change q/ t per time step is attributed to only E or to P .Corresponding to these assumptions, an identified decrease of specific humidity is attributed to the formation of precipitation.In cases of a positive increment of specific humidity, the method assumes moisture uptake from evaporation at ground level.In this case, the contribution of surface evaporation (f m ) in a time interval m to total specific humidity at the end of this time interval (q m ) is The formation of precipitation in a later time interval does not affect the f m calculated for the earlier time interval.In contrast to that, further moisture uptake in a later time interval n reduces the relative contribution of surface evaporation in the earlier time interval m to the moisture at the end of n.
In this case, the contribution from earlier moisture uptake to q n is recalculated according to In cases of moisture uptake above the ABL an attribution to surface evaporation is not directly evident.A maximum altitude for consideration of moisture uptake might therefore be appropriate.However, Aemisegger et al. (2014) relate moisture uptake at higher levels for trajectories starting from Rietholzbach in northern Switzerland to the outflow of shallow convection.As moisture in Karlsruhe and Rietholzbach originates from similar source regions, we also do not assume a maximum altitude for the consideration of surface evaporation.

Lagrangian isotope model
A Lagrangian isotope model was developed to serve as a benchmark for our δD measurements.Like the Lagrangian moisture source diagnostic, the model runs along kinematic 5-day back trajectories and attributes changes of specific humidity (q) to the formation of precipitation or moisture uptake from surface evaporation.Analogous to the moisture source diagnostic, the model does not apply a maximum altitude for moisture uptake, smoothes q along the trajectories with a 24 h rectangle function, and uses the same trajectory ensembles.From each trajectory ensemble, nine modeled δD values for Karlsruhe are obtained, which are combined to one average value by weighting the nine values with q at the arrival.

Dehydration
A decrease of specific humidity in a time interval indicates the formation of precipitation.Because of preferential fractionation of D into the liquid phase, the formation of precipitation results in a decreasing isotope concentration ratio  (Rayleigh and Ramsay, 1894).Based on a fractionation factor α D (T ), this Rayleigh model calculates changes of R D for infinitesimal changes of specific humidity q: (4) Specific humidity and air temperatures along the back trajectories were derived from the identical GDAS data set used for the calculation of the back trajectories.Under equilibrium conditions α D only depends on the temperature (T ) of an air parcel and increases from about 1.082 to 1.240 between +20 and −40 • C. For T >=0 • C we use a parameterization of α D over liquid water (Horita and Wesolowski, 1994).For T < 0 • C, we assume enhanced fractionation over ice and use the parameterization of Jancso and Van Hook (1974).A reduction of the fractionation factor α D (T ) in the case of supersaturation in ice clouds (Merlivat and Jouzel, 1979;Jouzel and Merlivat, 1984) by a factor on the order of 0.964 to 1 between −40 and 0 • C was considered following Jouzel and Merlivat (1984), with a supersaturation parameter λ of 0.004 according to Risi et al. (2010).

Moistening
In the case that specific humidity increases in a time interval [t 1 , t 2 ], we assume a permeable air parcel which takes up moisture q by turbulent mixing.We attribute that moisture to evaporation at ground level with the isotope ratio R D,E , which was lofted via small-scale turbulence or convection to the trajectory level.To calculate the corresponding change of R D of moisture in the tracked air parcel, we apply the following mixing equations: q(t 2 ) = q(t 1 ) + q, (5) If the entrained moisture from evaporation was transported via small-scale turbulence to the trajectory altitude, mixing with air from below the trajectory level is likely.However, applying the above equations, changes of R D and q in consequence of mixing with low-level air masses are ignored.We therefore implicitly assume that the air masses below the trajectory level experienced a similar transport and precipitation history as the tracked air parcel.R D and q of the tracked air parcel are not affected by mixing with that air from below, only the by freshly evaporated moisture.
Depending on the type of ground and skin temperature (T skin ), we calculate R D,E assuming evaporation from the ocean (R D,E_ocean ), continental evapotranspiration (R D,E_ET ), evaporation of melted snow (R D,E_snowevap.), or sublimation of snow (R D,E_snowsubl.).
For evaporation from the ocean, we assume vapor pressure fractionation over liquid water with α D (T skin ) according to the parameterization of Horita and Wesolowski (1994) and T skin at the trajectory position.To account for additional kinetic fractionation during evaporation on the order of α D,kin = 1.002-1.007(Merlivat and Jouzel, 1979;Pfahl and Wernli, 2009) due to different coefficients of diffusion of the different water isotopologues, we increased the α D by the factor α D,kin = 1.005.Please note that α D,kin is much smaller than the vapor pressure fractionation factor α D .Therefore, considering a dependence of α D,kin on environmental conditions is less important for R D than for other frequently used isotope ratios such as [H 18  2 O] / [H 16 2 O].Given the above assumptions, the isotope ratio of evaporated moisture R D,E_ocean only depends on the surface temperature and the isotope ratio of sea surface water (R D,ocean ) at the trajectory position: We derived R D,E_ocean from observations of δD and δ 18 O of ocean surface water collected in the Global Seawater Oxygen-18 Database (Schmidt et al., 1999).Since little data with the δD of ocean surface water exist, we calculated a median δD/δ 18 O ratio of 6.56 from Fröhlich et al. (1988), Duplessy (1970), Delaygue et al. (2001), Gat et al. (1996), Ostlund et al. (1987), Aharon and Chappell (1986), Yobbi (1992), and Weiss et al. (1979) and used this ratio to calculate δD from δ 18 O data.The δ 18 O of ocean surface water along the back trajectories was derived from the spatial 1×1 • interpolation (V1.1) of the Global Seawater Oxygen-18 Database by LeGrande and Schmidt (2006) (Fig. 1a), which we linearly interpolated to the locations along the trajectories.To derive skin temperatures representative of conditions during maximum evaporation (T skin ), we weighted skin temperatures along the trajectories (T skin,unweighted ) with positive surface latent heat flux in time intervals of ±2 h (Fig. 2).If there were less than 12 trajectory points (time resolution of 1 h) with significant latent heat flux above 2 W m −2 in an interval, it was extended for ±1 h until it contained 12 data points.T skin,unweighted and accumulated surface latent heat fluxes along the back trajectories were derived from a reduced GDAS data set with the same horizontal resolution of 1 × 1 • as it was used for the calculation of the back trajectories but with data only every 6 h.The data were interpolated linearly in space and time to the locations along the trajectories.The accumulated surface latent heat flux from the GDAS was divided by 6 to account for the hourly resolution of the trajectories.
Over the continent, evapotranspiration consists of evaporation from the bare soil, transpiration of plants, and evaporation from canopy interception.As a first simplification we ignore canopy interception; i.e., we consider condensation with subsequent complete re-evaporation as a neutral pro- cess with respect to q and R D .Moisture from the two other sources strongly differs in isotopic composition.
Evaporation from the bare soil is accompanied by isotope fractionation (Zhang et al., 2010).To calculate the isotope ratio of moisture evaporated from the bare soil (R D,E_soilevap.), we assume vapor pressure fractionation over liquid with α D (T skin ) according to the parameterization of Horita and Wesolowski (1994).To account for kinetic fractionation dur- ing evaporation from the soil on the order of α D,kin = 1.017-1.025(Mathieu and Bariac, 1996), we increased α D by the factor 1.021.We further assume the isotope ratios of soil water to be the same as the isotope ratios of precipitation (R D,prec.): Observations of R D,prec.have been collected in the Global Network of Isotopes in Precipitation (GNIP) (Araguas et al., 1996) of the IAEA since the 1960s.We used climatological monthly means of R D,prec. of the Regional Cluster-based Water Isotope Prediction (RCWIP) (Terzer et al., 2013), which provides a spatial interpolation of the GNIP data.RCWIP data are available with a horizontal resolution of 0.17×0.17• (Fig. 1b, c) and were linearly interpolated to the locations along the trajectories.Because soil water in central Europe is generally frequently recharged by precipitation, we ignore systematic enrichment of HD 16 O in the uppermost soil layer caused by continuous fractionating evaporation from the soil.For instance, measurements of precipitation amount at the measurement site in Karlsruhe indicate recharge of soil water around the site by more than 1 mm precipitation per day on average every 2.9 days.The assumption regarding enrichment is supported by findings of Risi et al. (2016), who observed insignificant systematic deviations of isotope ratios of water within the upper 15 cm of the soil from isotope ratios in precipitation at several sites in France, Germany, and the Czech Republic.
In contrast to bare soil evaporation, plants take up soil water from the soil and transpire that water completely, and therefore there is no fractionation into the atmosphere.Still some fractionation is possible on short timescales due to asynchronous accumulation and the release of HD 16 O and  H 16  2 O in leaves (Zhang et al., 2010).However, this process averages out over a day (Harwood et al., 1999;Farquhar et al., 2007) and is therefore ignored for modeling isotope ratios along the 5-day back trajectories.To calculate the isotope ratio of moisture originating from plant transpiration, we assume The isotope ratio of total evapotranspiration (R D,E_ET ) depends on the fraction of plant transpiration (FT) on total evapotranspiration: FT varies with region and on seasonal, synoptic, and diurnal timescales.For modeling we ignore these variations and use a constant fraction of transpiration.Based on Choudhury et al. (1998), Lawrence et al. (2007), and Aemisegger et al. (2014) we assume an average FT in Europe of 0.7.
Whenever we observe moisture uptake at continental skin temperatures below 0 • C, we ignore transpiration of plants (FT = 0) and attribute the moisture to the evaporation of melted snow or ice.In this case, we again assume equilibrium fractionation over liquid.We further assume the isotope ratio of snow to be the same as the climatological monthly means of the RCWIP: In Sect. 5 we investigate a possible role of snow sublimation.In this case we define a skin temperature T subl,max .Below that temperature we assume complete layer-by-layer sublimation of snow without isotope fractionation:

Initialization
To initialize R D of air masses originating from altitudes below 2 km above ground level, we assume isotope ratios (R D,ini ) in a convectively well-mixed ABL, where water vapor and ocean surface water or soil water are in isotopic equilibrium: For initialization at surface temperatures above 0 • C, we calculate fractionation factors according to the parameterization of Horita and Wesolowski (1994).At skin temperatures In a low-pressure system near Iceland the tracked air parcel ascended to an altitude of 3200 m.During the last 3 days of transport to Karlsruhe the air parcel was sinking to the sampling altitude.(b) After initialization in the MBL, specific humidity q (light blue and green colored line) and δD values (thick black line) of the tracked air parcel were slightly decreasing due to the formation of precipitation (dashed light blue lines) within the first day.More pronounced formation of precipitation, in consequence of lofting in a low-pressure system near Iceland, resulted in a second decrease of q and the modeled δD value dropped accordingly.Due to moisture uptake (green lines) related to a descent of the air parcel in the subsequent days, q and the δD value increased until the air parcel reached Karlsruhe.Thin black curves illustrate the modeled δD for different initializations of δD (Sect.2.3.3).The dependence on the initialization decreases with the amount of moisture uptake along the trajectories and is only low in Karlsruhe.
below 0 • C we assume fractionation over ice and apply the parameterization of Jancso and Van Hook (1974).
For air masses originating from an altitude (h) higher than 2 km above ground level we assume a linear decrease of R D,ini from boundary layer ratios at an altitude of 2 km to R D,10 km of 0.45 (Hanisco et al., 2007;Sayres et al., 2010) at an altitude of 10 km: The dependence of modeled isotope ratios from initialization decreases with moisture uptake along the back trajectories (Fig. 3).The uncertainty of R D,ini is especially strong at high altitudes, where air masses are strongly dehydrated and have a long history of isotope fractionation.Because air masses from high altitudes take up a lot of humidity during descent and transport to Karlsruhe, the uncertainty of R D from initialization is strongly reduced in Karlsruhe.Back trajectories corresponding to smaller moisture uptake typically originate from the ABL.In Karlsruhe R D of these back trajectories depends more strongly on R D,ini , which is much better defined within the ABL.Considering typical atmospheric moisture residence times in the range of 4-8 days (Läderach and Sodemann, 2016;Trenberth, 1998), using 10-day back trajectories instead of 5-day back trajectories should almost eliminate model uncertainty of the δD at Karlsruhe from uncertainty of R D,ini .However, such long trajectories easily cover distances of 10 000 km, which makes the modeled δD sensitive to potentially very different conditions in distant regions.This, in turn, makes uncertainty assessment of the modeled δD more complex.Using 5-day back trajectories is therefore a trade-off between a reasonably small sensitivity of the modeled δD at Karlsruhe on R D,ini and concentrating the analysis to the North Atlantic and Eurasia.Whenever it is necessary for the interpretation of our results, we assess uncertainty of the modeled δD from the model initialization by changing R D,ini in different model runs.

Isotope water vapor measurements
The concentrations of H 16  2 O and HD 16 O in low-level water vapor were measured for 17 months on a research campus 12 km north of Karlsruhe in southwestern Germany (49.10 • N, 8.44 • E; 110.4 m a.s.l.).
For the continuous measurements we used a Picarro water isotopologue analyzer L2120-i, which analyzed the ambient water vapor with a sampling rate of 0.6 Hz.The measurement technique is based on cavity ring-down spectroscopy, where the beam of a tunable diode laser is directed through a cavity, filled with the air to be analyzed.Based on the ring-down time of the laser light intensity, absorption spectra are measured between 7183.5 and 7184 cm −1 .A characterization of two similar analyzers (L1115-i and L2130-i) can be found in Aemisegger et al. (2012).Please note that the isotopologue analyzer also measures concentrations of the water isotope H 18 2 O.As the isotope concentration ratio  tube with a diameter of 6.4 mm.To reduce wall effects, we permanently flushed the inlet line with 30 standard L min −1 of air and used tubing made of electropolished stainless steel.To avoid condensation, the wall temperature of the 5 m of tubing inside the building was regulated to 22 • C. Saturation humidity corresponding to this temperature is above the analyzer limit of 14.9 g kg −1 .On 5 days in August 2012, humidity slightly exceeded that analyzer limit.Corresponding measurements were removed from the time series.
We report isotope measurements in the δ-notation, which normalizes isotope ratios to a standard scale, defined by the Vienna standard mean ocean water (VSMOW: δD = 0 ‰) and standard light Antarctic precipitation (SLAP: δD = −428.0‰) (IAEA, 2009).For the automated calibration of the analyzer we applied a Picarro standard delivery module (A0101), which allows the alternating injection of two different water standards into a Picarro vaporizer (A0211).In this vaporizer the liquid standards immediately evaporate in a constant flow of dry synthetic air (140 • C, 0.3 standard L min −1 dry air flow with 1.2 mg kg −1 residual humidity).A two-point calibration was done for 2 h every 10 h at δD of −62.1 and −142.2 ‰.
Instrumental drift during the 17 months was below 3 ‰ (Fig. 4a).The total accuracy of our δD measurements due to uncertainty of calibrations and instrumental drift between two calibrations was 0.98 ‰.The 0.6 Hz precision of the www.atmos-chem-phys.net/17/1207/2017/measurements was below 1 ‰ and can be ignored for 10 min averages shown in this paper.
Based on the two-point calibrations, we applied linear stretching to the measurements.Sixty-three percent of our observations are within the range of isotope ratios covered by the two standards.To approve linearity of the applied correction for isotope ratios below that range, we performed repeated calibrations with a third standard (δD = −245.3‰) in the 2 years subsequent to the campaign.We found additional uncertainty of δD at −245.3 ‰ due to a slight nonlinearity of the applied correction to be smaller than 0.3 ‰, which is in agreement with the more detailed characterization of Aemisegger et al. (2012).
To identify a potential humidity dependence of the isotope ratio measurements, we generated three humidity levels between 1.8 and 13.7 g kg −1 during each calibration.We found the humidity dependence of the instrument to be smaller than the uncertainty of individual calibrations.Therefore, we only applied the average humidity dependence found using all calibrations of −0.021 ‰ (g kg −1 ) −1 to the data set (Fig. 4b).

Meteorological data at the measurement site
Observations of specific humidity at the measurement site were derived from the Picarro isotope analyzer.For calibration of the Picarro humidity measurements we used observations of a VPT6 Thygan dew point mirror hygrometer (Meteolabor, Switzerland), which was mounted on a meteorological tower 30 m above ground level 900 m in the westsouthwest.Since the topography at the measurement site is flat for some kilometers in all directions, we assume the tower observations to be representative for the measurement site.The dew point hygrometer performed a measurement of 1 min every 10 min and has an uncertainty of ±0.1 K. Tenminute averages of specific humidity derived from the Picarro Analyzer and observations of the dew point hygrometer show a correlation coefficient R of 0.9913.For the calibration of the Picarro humidity observations we applied the mean linear regression between hygrometer data and Picarro measurements of 1.13, which varied by 1 % between the first and second half of the measurement period.
The amount of precipitation was measured at the meteorological tower at ground level with a time resolution of 10 min.

Analysis of seasonal and synoptic variations
In this section, we present the measurements from Karlsruhe, covering the time period from January 2012 to May 2013.For this time period, we identify specific circulation regimes related to cold snaps in Karlsruhe.Subsequent to this, we examine the capability of the Lagrangian isotope model of reproducing corresponding variations of δD.

Continental temperatures and zonal circulation
According to 5-day back trajectories, 85 % of the sampled low-level air masses at Karlsruhe originated from altitudes below 2 km above ground level.Consequently, most of the air masses were exposed to continuous moisture uptake from surface evaporation during transport to Karlsruhe.For identification of major source regions of the water vapor in Karlsruhe we applied the Lagrangian moisture source diagnostic described in Sect.2.2.Based on this analysis, we found the North Atlantic to be the most important moisture source (Fig. 5), from where westerlies transported the tracked air masses to the measurement site.In addition to the predominantly westerly moisture transport, inversions of zonal circulation in winter occasionally led to easterly moisture transport.Because of the finite length of the 5-day back trajectories, the total identified humidity is lower than 100 %.When humidity is smoothed along the back trajectories for 24 h, the total identified humidity accounts for 47 %.If smoothing specific humidity for 12 h, diurnal and sub-diurnal variations in humidity are interpreted as the formation of precipitation and moisture uptake, which increases the total identified humidity to 63 %.Both numbers are in reasonable agreement with more extended studies implying global atmospheric moisture residence times in the range of 4-8 days (Läderach and Sodemann, 2016;Trenberth, 1998).
Air temperatures, specific humidity, and δD in Karlsruhe followed similar seasonal patterns (Fig. 6).In winter (December, January, February), air temperatures 30 m above ground level (T ) were on average 2 • C. Towards summer (June, July, August), T increased to on average 20 • C. Higher T in summer corresponds to a higher saturation vapor pressure.That allows the transport of marine air to Karlsruhe with less condensation in summer than in winter.Consequently specific humidity (q) in Karlsruhe rose from 6 g kg −1 (DJF) to 14.8 g kg −1 (JJA).δD changed from −162 ‰ (DJF) to −109 ‰ (JJA) and thereby consistently with T and q implies a lower degree of condensation and rainout in summer than in winter.
The gray color in Fig. 6 identifies circulation regimes with easterly moisture transport (25 % of data).Such regimes predominantly occurred in winter and resulted in the transport of continental air masses to Karlsruhe.The corresponding air masses usually were marked by especially low T and q, which led to pronounced cold snaps in Karlsruhe.Consistent with the low T and q, the air masses during cold snaps showed an especially low δD value.
Both findings -the seasonality of δD and the especially low δD in cold, continental air masses from the east -are in good agreement with the well-known "continental effect".That effect describes a decrease of δD in precipitation over continents with distance to the coast (Fig. 1b), caused by the relation between δD and degree of rainout.Since the δD of rain depends on the δD of the water vapor it is formed from, it is reasonable to find a similar continental effect imprinted to water vapor as well.

Comparison of measured and modeled δD
The Lagrangian isotope model described in Sect.2.3 is able to reproduce the observed slow seasonal variation of δD, as well as the strong and relatively fast variations of δD due to circulation regimes with predominantly easterly moisture transport in winter (Fig. 7a, b).The mean difference between modeled and observed δD ( δD) is −4.1 ‰.The correlation coefficient R of modeled and observed δD is 0.82.Thereby, the correlation calculated for different seasons strongly differs from the overall correlation.For summer R is only −0.03 due to the small variability of δD in this season.For winter R is 0.87.
Figure 7c shows the scatter plot of measured and modeled δD.Furthermore, the figure illustrates the impact of the formation of precipitation and surface evaporation on the modeled δD.If the formation of precipitation is ignored in the model, the overall correlation with the observations is still 0.79.The main reason for this is the relation of observations with low δD values to easterly moisture transport.Corresponding back trajectories are initialized with relatively low δD values according to GNIP observations in the respective continental source regions.This means that we only ignore the formation of precipitation in the 5 days covered by the back trajectories, but we implicitly consider the formation of precipitation which determined the isotope ratios in precipitation in the moisture source regions.However, due to surface evaporation with relatively high δD values, the δD in this scenario is +23.9 ‰.If one considers the formation of precipitation but ignores the surface evaporation, the overall R is reduced to 0.62.The corresponding δD is −34.8 ‰.So consideration of both processes, surface evaporation and the formation of precipitation, is essential for reproducing the observed δD.
With respect to surface evaporation the air masses from the west and from the east contain very different types of information.The air masses from the west are exposed to moisture uptake from the ocean and to continental evapotranspiration at relatively warm temperatures.Therefore they contain information about isotope fractionation during evaporation from the ocean, evaporation from warm land surfaces, and plant transpiration.This makes the modeled δD of westerly air masses especially sensitive to simplifying model assumptions regarding the δD of moisture from continental evapo-transpiration at warm skin temperatures.Increasing for instance the fraction of plant transpiration on total evapotranspiration from 0.7 to 0.8 in the model increases the mean modeled δD from summer by 2.0 ‰.Assuming δD values of soil water which are systematically increased by +10 ‰ increases the mean modeled δD from summer by +2.4 ‰.In contrast to the westerly air masses, the easterly air masses during cold snaps are sensitive to isotope fractionation during surface evaporation at temperatures where snow and meltwater exist.For these easterly air masses assumptions regarding evapotranspiration at warm skin temperatures only have a small impact.
An interesting feature in Fig. 7a are spikes of low δD values (blue), which are not reproduced by the model (red).We frequently observed such spikes from spring to autumn.Potential processes causing the spikes are the evaporation of rain below the cloud base or isotope exchange between falling raindrops and the low-level water vapor.An impact of these sub-cloud processes on the δD of low-level water vapor is demonstrated for individual weather fronts by Wen et al. (2008) and Aemisegger et al. (2015) and is complementarily supported by observations of isotope ratios in precipitation (Friedman et al., 1962;Stewart, 1975;Gedzelman and Arnold, 1994).As isotope processes below clouds are not represented in the Lagrangian isotope model, the observations related to sub-cloud processes cannot be further investigated by means of this model.However, a relevant role of sub-cloud processes for our observations of δD spikes is supported by the strongly increased probability of precipitation during the spikes.For respective observations with δD values smaller than expected from the standard deviation between observed and modeled δD (23.5 ‰), the probability to observe precipitation in Karlsruhe within ±3 h was 44 %, whereas for the other observations this probability was only 24 %.For the investigation of δD during cold snaps in winter sub-cloud processes can be ignored, because the interaction between falling precipitation and water vapor is strongly suppressed in cases of solid precipitation.

δD during cold snaps
The good agreement of modeled and measured δD in winter (R = 0.87) underlines the strong potential of the Lagrangian isotope model for analyzing isotope processes in the remote moisture source regions during cold snaps.In this section, we select respective observations and investigate isotope fractionation during sublimation or evaporation at ground level.

Sublimation of snow or snowmelt evaporation?
Evaporation below 0 • C was historically often considered as non-fractionating layer-by-layer sublimation of snow and ice (Ambach et al., 1968;Dansgaard, 1973;Friedman et al., 1991) because of a low coefficient of self-diffusion of wa- ter molecules in ice.However, this assumption ignores snowmelt and fractionation during evaporation from the liquid phase, as implied by Gurney and Lawrence (2004), Lechler andNiemi (2011), andNoone et al. (2013).Both assumptions result in very different δD of moisture from evaporation.In the case of non-fractionating sublimation, this δD equals the δD of the snow.In the case of fractionating evaporation, the δD value of moisture from evaporation is about 90 ‰ lower.Our trajectory model provides an opportunity to test both formulations and to assess which one is more reasonable for central Europe.
According to the Lagrangian moisture source diagnostic, only 5 % of the humidity analyzed in Karlsruhe from January 2012 to May 2013 originated from surface evaporation at locations with a GDAS-based snow depth greater than 0.5 cm.For this reason, we used data fulfilling three selection criteria for the further investigation of isotope fractionation during evaporation from snow-covered surfaces.First, we excluded observations which are strongly affected by evapotranspiration at skin temperatures (T skin ) above the freezing point.For this purpose, we identified moisture uptake at T skin > 0 • C by means of the moisture source diagnostic and excluded air masses with a respective contribution above 2 %.For air masses meeting the first criterion the median contribution of moisture evaporated at T skin < 0 • C is 28 %.The air masses with the smallest moisture uptake within the last 5 days are least sensitive to the δD of moisture from surface evaporation and therefore do not allows us to robustly evaluate the model description of isotope fractionation during the sublimation of snow and the evaporation of meltwater.For a meaningful interpretation, we therefore only used the half of data with a contribution from surface evaporation at T skin < 0 higher than 28 %.Only 2 % of the corresponding air masses originated from altitudes higher than 2000 m above ground level.Since uncertainty of model initialization is especially high for these air masses, we finally excluded the 2 % of air masses originating from high altitudes.
Some 174 of the 3-hourly modeled data points meet the three selection criteria.They belong to 38 different days and were used for further interpretation.Respective air masses mainly originated from the east (Fig. 8a).The mean fraction of moisture identified for these air masses by the moisture source diagnostic is 48 %.When smoothing humidity along the back trajectories for 12 h instead of 24 h, the mean identified fraction is 68 %.The GDAS data indicate the existence of snow on the ground at 96 % of locations along the selected back trajectories (Fig. 8b) with a median snow depth of 1.8 ± 0.9 cm (± gives the interquartile range).The median T skin at the trajectory positions 5 days back was −11.7 ± 5.1 • C.During transport to Karlsruhe the T skin rose on average by 6.0 ± 3.3 K.A decrease of relative humidity due to warming of the air masses was partially compensated by moisture uptake and a corresponding increase of specific humidity by on average 52 ± 35 %.
Assuming equilibrium fractionation during evaporation of meltwater at T skin < 0 • C in the reference run (M MW , Table 1), the model underestimates the selected δD values by on average δD = −18.6±1.5 ‰ (± gives the statistical uncertainty of the mean).Assuming non-fractionating sublimation at T skin < 0 • C in a further model run (M S , Table 1) results in δD values that are on average δD = +26.9± 1.6 ‰ above the observations.
Considering the relatively high δD values of moisture from non-fractionating sublimation and the about 90 ‰ lower δD values in the case of fractionating evaporation, the difference of mean δD between both model runs is qualitatively reasonable.To judge whether one run provides more realistic results, we tested if one of the two scenarios could be brought into agreement with the observations when considering uncertainty of model assumptions.
As the most important source of uncertainty for modeling the δD in Karlsruhe during cold snaps we consider variability of the δD of surface layer snow in the moisture source regions that is not fully captured by the model.This variability may systematically affect (1) the assumed δD of moisture from surface evaporation as well as (2) the δD at the model initialization.In addition, (3) the amount of identified moisture uptake needs to be accurate to reliably simulate the impact of surface evaporation on the δD of water vapor.
In the following, we estimate uncertainty of respective model assumptions.Subsequent to this, we vary the model Table 1.Average differences between modeled and measured δD of the 174 selected data points ( δD), data points of group "cold" ( δD cold ), and data points of group "warm" ( δD warm ) from different model runs (M).± states the statistical uncertainty of the averages (root mean square error divided by the square root of the number of observations).Values of particular interest are printed in bold type.
1.The δD of moisture from sublimation or evaporation of meltwater depends on the δD of the surface layer snow, which we assume to be equal to the δD of precipitation.
The δD of precipitation in an individual year may systematically differ from the used climatological monthly means of the RCWIP.To assess typical interannual variability of δD of precipitation in winter, we used data from 134 European and Russian GNIP stations from the midlatitudes between 8.4 and 50 • E with observations from at least 3 years.For each of the stations we calculated the mean δD in the different winters (November, December, January, and February) and the standard deviation of the winter averages.The mean of standard deviations of the different stations was 11.5 ‰, which we assume to reflect the mean interannual variability of δD in precipitation in the moisture source regions related to cold snaps.
Because there is a general relation between surface air temperatures and δD of precipitation in central Europe (Schoch-Fischer et al., 1983;Jacob and Sonntag, 1991), winter months from years with especially high air temperatures and a potentially strong contribution from liquid precipitation are related to relatively high δD values.Since we want to estimate the δD of snow, data from winters with especially low temperatures and a strong contribution from solid precipitation with low δD values are likely to be most representative.During these winters the δD of precipitation is probably closest to δD of the RCWIP minus the 11.5 ‰.
In contrast to this, post-depositional fractionation processes may increase the δD values of surface layer snow on the order of +10 to +20 ‰ during periods with small accumulation rates (Gurney and Lawrence, 2004;Moser and Stichler, 1974), suggesting a scenario where the δD values of surface layer snow are higher than the δD from the RCWIP.Please note that this scenario is not likely for a seasonal snowpack during melt season since ablation may uncover old snow from colder winter months, causing changes of the δD of surface layer snow of −50 ‰ (Dahlke and Lyon, 2013).Given the relatively small snow depths in the investigated moisture source region (Fig. 8b), we ignore the uncovering of older snow layers and only consider a potential postdepositional increase of the δD of surface layer snow.
To test whether the too high modeled δD values from the scenario of sublimation (M S ) can be significantly reduced when considering systematic uncertainty of the δD of surface layer snow regarding interannual variability of the δD of snowfall, we performed one model run (M S-,snow ), in which we shifted the δD of snow by −11.5 ‰.To test whether the too low modeled δD values from the scenario of evaporation of meltwater (M MW ) can be increased by considering a postdepositional increase of the δD of surface layer snow, we performed one further model run (M MW+,snow ), in which we shifted the δD of surface layer snow by the same absolute value of +11.5 ‰.The mean difference between modeled and observed δD ( δD) of the different model runs is listed in Table 1.
2. During cold snaps in Karlsruhe, on average 48 % of humidity could be attributed to moisture uptake along the 5-day back trajectories.The other side of this argument is that the δD of 52 % of humidity in Karlsruhe is determined by the initialization of isotope ratios.
For initialization, we assume δD in a well-mixed ABL and isotopic equilibrium between water vapor and the climatological monthly δD of precipitation.This assumption is in agreement with one of the rare extended, simultaneous time series of δD in water vapor and precipitation, conducted 45 km NNE from our site (Jacob and Sonntag, 1991).This study shows monthly averages of δD in precipitation and water vapor at ground level for the years 1981-1988.The average deviation of the δD of water vapor to isotopic equilibrium with precipitation in November, December, January, and February was +4.1 ± 7.7 ‰ (± states the standard deviation calculated from the winter averages of the different years).
To test how much the values of δD for M S and M MW can be reduced by considering the uncertainty of δD at the model initialization, we performed two further model runs.For the scenario of sublimation we performed a model run in which we shifted δD at the initialization for −3.6 ‰ (M S-,ini ).For the scenario of evaporation of meltwater we shifted δD at the initialization for +11.8 ‰ (M MW+,ini ).
3. Potentially artificial diurnal variations of q along the back trajectories could result in an overestimation of the formation of precipitation and moisture uptake.To avoid such an overestimation, diurnal variations of q were suppressed by smoothing q along the back trajectories with a 24 h broad rectangle function.Arbitrarily choosing a width of 24 h may smooth out real subdiurnal and diurnal variations of q and thereby may result in an underestimation of the amount of moisture uptake.To assess the potential impact of the smoothing on the modeled δD, we changed the width of the applied rectangle kernel to 12 h in M S-,upt.12h and to 36 h in M MW+,upt.36 h .
To finally assess the minimum possible values of δD in the case of superposition of the three sources of uncertainty discussed above, we combined the assumptions of M S-,snow , M S-,ini , and M S-,upt.12h in the model run M S---and M MW+,snow , M MW+,ini , and M MW+,upt.36 h in the model run M MW+++ .
Table 1 summarizes δD for the different model runs.None of the model runs considering only one source of uncertainty is able to reduce δD for the scenarios of sublimation or evaporation of meltwater to values close to 0. Even for M S---, which simultaneously assumes all the uncertainties of model assumptions in the scenario of sublimation, δD is +16.0 ‰.The discussed uncertainty terms are therefore not able to bring model and observations into agreement with each other when only considering non-fractionating sublimation.This implies that fractionating evaporation of meltwater played a significant role during our observations.
For M MW+++ , which simultaneously assumes all the uncertainties of model assumptions in the scenario of evaporation of meltwater, the absolute value of δD is reduced to 0.9 ‰.Considering the statistical uncertainty of δD of 1.6 ‰, the average modeled and measured δD of the 174 selected air masses may therefore be brought into rough agreement with each other when assuming fractionating evaporation of meltwater.However, this requires superposition of the different uncertainty terms.
In order to refine this result, we split the selected observations into two groups of equal size according to the predominant skin temperature during moisture uptake (T skin,predom.).For this purpose, we weighted skin temperatures along the individual ensembles of back trajectories with moisture uptake identified by the Lagrangian moisture source diagnostic.The median T skin,predom. of the 174 selected trajectory ensembles is −6.92 • C. We attributed data points to a group "cold" if T skin,predom. of the respective trajectory ensemble is below −6.92 • C. For T skin,predom.above −6.92• C we attributed data points to a group "warm".Please note that also points of group "warm" have a T skin,predom.below 0 • C according to our selection criteria.Due to interannual variability of T skin,predom., data of the two groups are not randomly distributed in time.Seventy-seven percent of group "cold" corresponds to an especially pronounced cold snap in February/March 2012, whereas 85 % of group "warm" belongs to cold snaps between October 2012 and February 2013.
Figure 9 shows two-dimensional probability distributions of the selected modeled and measured δD.Blue denotes data from group "cold" and red denotes data from group "warm".Under the assumption of non-fractionating sublimation (M S ), modeled δD values of both groups are significantly higher than the observed values (Fig. 9a).The overestimation of modeled δD values is especially strong in the regime with higher T skin , where snowmelt and fractionating evaporation would be expected.Figure 9b shows the respective probability distributions under the assump-tion of snowmelt and fractionating evaporation of meltwater (M MW ).Under this assumption, the modeled δD of group "warm" is close to the observations.However, modeled δD values of group "cold" are now far too low.
Table 1 lists the mean differences between modeled and measured δD for the two groups ( δD cold , δD warm ).As δD of the individual groups may deviate more from the observations than the mean δD of all 174 selected air masses, analyzing δD cold and δD warm allows drawing less ambiguous conclusions than analysis of δD.For M S,---the value of δD warm (18.2 ‰) is larger than the value of δD of all selected data (16.0 ‰), which underlines the importance of fractionating evaporation for reproducing the observations.For M MW,+++ the absolute value of δD cold (10.5 ‰) is larger than the absolute value of δD (0.9 ‰).So even M MW,+++ , in which we simultaneously assumed all the uncertainties of model assumptions in the scenario of fractionating evaporation of meltwater, does not allow reproducing the observations of group "cold".This, in turn, implies significant non-fractionating sublimation during our observations.
Comparison of δD observations with δD of the Lagrangian isotope model therefore implies a relevant role of both types of isotope fractionation in central Europe.

Temperature-dependent types of fractionation
To simultaneously bring into agreement modeled and observed δD of data corresponding to group "cold" and group "warm", we suggest the existence of two regimes of T skin with predominant non-fractionating sublimation in the colder regime and predominant fractionating evaporation of meltwater in the warmer regime.
For the characterization of these two regimes we assume a maximum temperature for non-fractionating sublimation in the model: T subl,max .For T skin < T subl,max , we assume nonfractionating sublimation.In the case of T subl,max < T skin < 0 • C, we assume equilibrium fractionation during the evaporation of meltwater.
To assess T subl,max for optimal agreement between modeled and observed δD, we performed 16 model runs with a different T subl,max in each run (−15 to 0 • C in steps of 1 K).We refer to these runs as M S_MW,T subl,max .The thick black line in Fig. 10 shows the mean differences between modeled and observed δD of the 174 selected data points ( δD) from the 16 different M S_MW,T subl,max .Modeled δD values are highest in M S_MW,0 • C , which assumes non-fractionating sublimation of snow for all T skin < 0   The thin dashed black line corresponds to the model configuration M S_MW,+++,T subl,max , which does not allow reproducing the observations in group "warm" and is therefore not considered in the gray shaded area of uncertainty.In magenta, optimal T subl,max for best agreement of model and observation (dot) and uncertainty of the optimal T subl,max (error bar) due to uncertainty of model assumptions are shown.
The agreement between observed and modeled δD from the M S_MW,T subl,max is best for a T subl,max of −7.7 • C. For this T subl,max the mean δD of all selected data points is 0 and also δD of group "cold" as well as the δD of group "warm" is approximately reproduced by the model.Figure 9c shows the respective two-dimensional probability distributions.
The statistical uncertainty of this optimal T subl,max due to scatter between modeled and observed δD is 0.7 • C. Further uncertainty is introduced to the determined optimal T subl,max by the assumptions of the Lagrangian isotope model, which can systematically change the mean modeled δD and the optimal T subl,max .To assess this uncertainty of the optimal T subl,max , we performed 16 • 8 = 128 additional model runs (Fig. 10, thin lines) corresponding to 16 different T subl,max from −15 to 0 • C and 8 different model configurations with the same assumptions about a changed δD of snow, a changed δD at the model initialization, a different amount of moisture uptake, and superposition of the three effects as in the above uncertainty assessment for the M S and the M MW .Analogous to the uncertainty assessment for the M S and the M MW , we refer to the model runs as M S_MW,+,snow,T subl,max , M S_MW,+,ini,T subl,max , M S_MW,+,upt.36 h,T subl,max , M S_MW,+++,T subl,max , M S_MW,-,snow,T subl,max , M S_MW,-,ini,T subl,max , M S_MW,-,upt.12 h,T subl,max , and M S_MW,---,T subl,max .The model runs with assumptions related to higher modeled δD values (dashed thin lines) result in a lower optimal T subl,max and model runs with assumptions related to lower modeled δD values (solid thin lines) result in a higher optimal T subl,max .
Thin black lines in Fig. 10 depict the maximum possible shift of the average modeled δD in Karlsruhe in the case of superposition of the examined model assumptions.The solid thin black line reflects a maximum unfavorable superposition of assumptions related to lower modeled δD values (M S_MW,---,T subl,max ).The M S_MW,---,T subl,max therefore allow to assess the upper bound of −3.6 • C of the T subl,max for optimal agreement between model and observation.The lower bound of the optimal T subl,max , which is derived from the M S_MW,+++,T subl,max , is −15 • C (thin dashed line).Here again, individual analysis of data from the groups "cold" and "warm" allows us to refine the result.For M MW,+++ δD warm is +8.6 ‰ (Table 1).Since M MW,+++ assumes fractionating evaporation of meltwater for all T skin < 0 • C, it marks the lower boundary of δD values from the M S_MW,+++,T subl,max .The corresponding set of assumptions does therefore not allow reproducing the observations in group "warm" even when assuming a very low T subl,max .Assuming maximum unfavorable superposition of the uncertainties of model assumptions in the M S_MW,+++,T subl,max is therefore too conservative for the uncertainty assessment of the optimal T subl,max .For this reason, we assess the lower bound of the optimal T subl,max by means of the M S_MW,+,snow,T subl,max , which only consider one uncertainty term and just allow the reproduction of the observed δD in group "warm".From the M S_MW,+,snow,T subl,max we derive a better confined lower bound of the optimal T subl,max of −11.2 • C.
So the uncertainty of model assumptions translates into an uncertainty range of T subl,max for optimal agreement of model and observation from −11.2 to −3.6 • C. Together with the statistical uncertainty of 0.7 • C, the total uncertainty range of T subl,max sums up to −11.9 to −2.9 • C.
In this paper, we investigated isotope fractionation during surface evaporation in snow-covered regions.For this purpose, we combined 17 months of measurements of δD in low-level water vapor in central Europe with a new Lagrangian isotope model.
By means of this approach, we identified two regimes of GDAS skin temperatures (T skin ) below the freezing point with significantly different deviation between modeled and observed δD.To resolve this difference, we suggest two regimes of T skin with different types of predominant isotope fractionation.Based on sensitivity tests with the Lagrangian isotope model, we found that the colder regime is described best by non-fractionating sublimation of snow.The warmer regime is described best by fractionating evaporation of meltwater.
We determined a T subl,max separating both regimes of T skin by optimizing the agreement between modeled and observed δD.For a T subl,max of −7.7 • C this agreement is best.Uncertainty related to assumptions of the isotope model corresponds to a range of uncertainty of T subl,max from −11.9 to −2.9 • C.
The finding of a cold temperature regime with a small impact of fractionation during sublimation at ground level does not contradict earlier studies on snow which indicate fractionating interaction between the surface layer snow and atmospheric water vapor, even in cases of temperatures far below the freezing point (Epstein et al., 1965;Stichler et al., 2001).These studies document systematic post-depositional increases of isotope ratios in the snowpack, which imply processes such as slight kinetic fractionation during sublimation in consequence of different coefficients of diffusion of the different water isotopes or fractionating vapor deposition.Given the uncertainty of the Lagrangian isotope model, these small effects would not be detectable by our approach.Nevertheless, these effects might result in significant postdepositional modifications of isotope ratios in the snowpack on timescales longer than the 5 days covered by the trajectories.
For GDAS skin temperatures between T subl,max and 0 • C our results imply significant fractionating evaporation of meltwater.The identification of a fractionating "meltwater regime" is consistent with earlier observations of isotope ratios in snow, which point to fractionating evaporation during ablation at temperatures close below the freezing point (Moser and Stichler, 1974;Gurney and Lawrence, 2004;Lechler and Niemi, 2011).Since snow samples give an integrated signal over long time periods, a detailed attribution of these observations to certain meteorological conditions is difficult.Complementary to the studies on snow, the method presented here allows the post-depositional isotope fractionation to be attributed to meteorological conditions with GDAS skin temperatures between T subl,max and 0 • C. The determined T subl,max refers to a GDAS skin temperature that was weighted with positive latent heat flux at ground level.For this reason T subl,max is representative for GDAS skin temperatures during the day, when evaporation is strongest.However, it should be kept in mind that due to the coarse resolution of 1 × 1 • of the GDAS data, much spatial variability of the skin temperature is smoothed out.The meaning of the T subl,max derived in this study is therefore an average temperature in a 1 × 1 • grid cell above which the evaporation of meltwater exceeds the amount of moisture from sublimation.A way to derive a more physical temperature separating the regimes of sublimation and meltwater evaporation could be using data with a higher horizontal resolution.This would not necessarily improve accuracy with respect to the back trajectories' positions but locations with enhanced skin temperatures and especially high amounts of surface evaporation would be more realistically represented, presumably resulting in higher evaporation weighted skin temperatures and a higher T subl,max .
In addition, higher horizontal resolution would allow to better account for spatial heterogeneity of the δD of surface layer snow in mountainous regions, for instance by also weighting the δD of the snow with the amount of surface evaporation.In this context please note that systematic uncertainty regarding the δD of surface layer snow turned out to be the main limitation for determining T subl,max with the approach presented here.For this reason regular analysis of the δD of surface layer snow samples, for instance at selected GNIP stations, would be a very desirable and efficient measure to reduce uncertainty of T subl,max .
Our results show that surface evaporation in the two identified regimes of skin temperature has a strong impact on the δD of low-level water vapor in central Europe.For isotope applications based on relations between δD and temperature, the consideration of the different types of isotope fractionation in both regimes is therefore of great interest.For instance, seasonal ablation in coastal regions of Greenland might systematically affect the relation between the δD of water vapor and dew point temperature over Greenland.Because such a seasonality may be different in climatologically different time periods, it may introduce uncertainty to temperature reconstructions from Greenland ice cores.
Furthermore, fractionating evaporation of meltwater will increase the δD value of the residual meltwater and in the case of recrystallization, the δD value of the snowpack.Ignoring fractionating evaporation therefore introduces uncertainty to a variety of isotope applications from reconstructions of paleotemperatures and paleotopography to studies on the formation of groundwater.The specification of a temperature regime with enhanced fractionation during evaporation may therefore help to identify, investigate, and reduce biases inherent to these applications.

Figure 1 .
Figure 1.(a) δD of ocean surface water derived from the interpolation (V1.1) of the Global Seawater Oxygen-18 Database by LeGrande and Schmidt (2006), assuming a constant factor of 6.56 between δD and δ 18 O in ocean surface water.(b) Climatological δD in precipitation in winter (DJF) from the Regionalized Clusterbased Water Isotope Prediction (RCWIP), which in turn is based on observations of the Global Network of Isotopes in Precipitation (GNIP).(c) Same as in (b) but for summer (JJA).

Figure 2 .
Figure 2. (a) Probability distributions of continental GDAS skin temperatures (T skin,unweighted ) (blue) and of skin temperatures weighted with the accumulated hourly latent heat flux at ground level within ±12 h (T skin ) (green).The occurrence of low temperatures is reduced as a consequence of the weighting.A peak around 0 • C becomes more clearly visible.(b) Illustration of the weighting algorithm for one exemplary back trajectory (arrival in Karlsruhe 4 May 2012, 21:00 UTC).

Figure 3 .
Figure 3. Illustration of the isotope modeling for one exemplary back trajectory (arrival in Karlsruhe on 18 March 2012, 00:00 UTC).(a) Altitude of the back trajectory (black) and terrain height (gray/blue).The isotope model was initialized at 80 • N in the marine boundary layer (MBL).In a low-pressure system near Iceland the tracked air parcel ascended to an altitude of 3200 m.During the last 3 days of transport to Karlsruhe the air parcel was sinking to the sampling altitude.(b) After initialization in the MBL, specific humidity q (light blue and green colored line) and δD values (thick black line) of the tracked air parcel were slightly decreasing due to the formation of precipitation (dashed light blue lines) within the first day.More pronounced formation of precipitation, in consequence of lofting in a low-pressure system near Iceland, resulted in a second decrease of q and the modeled δD value dropped accordingly.Due to moisture uptake (green lines) related to a descent of the air parcel in the subsequent days, q and the δD value increased until the air parcel reached Karlsruhe.Thin black curves illustrate the modeled δD for different initializations of δD (Sect.2.3.3).The dependence on the initialization decreases with the amount of moisture uptake along the trajectories and is only low in Karlsruhe.
more sensitive to kinetic fractionation than R D , it would be essential to more detailed consider kinetic fractionation during evaporation for modeling R18 O with the Lagrangian isotope model.Uncertainty of the modeled R18 O related to uncertainty of the kinetic fractionation factor does not allow deeper insight from analysis of R D and R18 O than from analysis of R D alone.Thus, we do not use the H 18 2 O measurements in this paper.The Picarro water isotopologue analyzer was located on the sixth floor of the Institute of Meteorology and Climate Research -Atmospheric Trace Gases and Remote Sensing of the Karlsruhe Institute of Technology.A downward-facing inlet funnel was installed 1 m above the edge of the roof, which corresponds to an altitude above ground level of 28 m.The connection to the inlet was established with a 6 m long

Figure 4 .
Figure 4. (a) Deviation of individual calibration measurements with the Picarro water isotopologue analyzer from isotope ratios of the liquid standards.Each point represents a calibration of 1 h.Dots: Standard 1; crosses: Standard 2. (b) Humidity dependence of the δD measurements.Long-term drift depicted in (a) was removed in (b)by subtracting the 5-week running average of calibrations.Red regression lines were calculated for both standards simultaneously by subtracting the mean difference between both standards.± gives the difference between the slopes calculated for the first and second half of the measurement period.

Figure 5 .
Figure5.Source regions of moisture (q) 30 m above ground level in Karlsruhe (black star) based on 5-day back trajectories for the time period January 2012 to May 2013.The color code indicates the contribution of different source regions to q in Karlsruhe in % per km 2 .Integration over the whole map gives the total identified humidity of 47 %.

Figure 6 .
Figure 6.Measurements in Karlsruhe from January 2012 to May 2013 30 m above ground level (10 min averages).Black: 5-day back trajectories originate from the west; gray: 5-day back trajectories originate from the east.(a) δD of water vapor.Gaps in the time series are caused by instrumental issues with analyzer and calibration device.(b) Air temperature (T ).(c) Specific humidity (q).

Figure 7 .
Figure 7. (a) Time series of measured (blue) and modeled (red) δD of water vapor in Karlsruhe.(b) Enlarged section of (a), which demonstrates the capability of the model of capturing the high variability of δD in winter.(c) Measured versus modeled δD.The 3hourly available modeled δD is compared to the temporally closest 10 min average of the δD observations.Gray: reference run, surface evaporation (E), and the formation of precipitation (P ) are considered (R = 0.82, δD = −4.1 ‰); green: only E is considered and P is ignored (R = 0.79, δD = +23.9‰); magenta: only P is considered and E is ignored (R = 0.62, δD = −34.8‰); black: 1 : 1 line.

Figure 8 .
Figure 8.(a) Lagrangian source region analysis of low-level water vapor in Karlsruhe for the observed cold snaps.The mean identified fraction of moisture along the 5-day back trajectories is 48 %.(b) Mean snow depth during the cold snaps based on GDAS data.
• C and is therefore identical with M S .Out of the 16 model runs M S_MW,−15 • C gives the lowest δD values.The δD values from M S_MW,−15 • C are close to δD from M MW , as most T skin along the back trajectories were above −15 • C, which means that almost no sublimation below T subl,max is considered.M S_MW,−29 • C gives the same results as M MW , because there was no moisture uptake identified for T skin < −29 • C.

Figure 9 .
Figure 9. Two-dimensional probability distributions of measured and modeled δD of low-level water vapor in Karlsruhe for selected cold snap events.Red: group "warm"; blue: group "cold".Probabilities were calculated for a 20 ‰ × 20 ‰ δD grid, smoothed with a 20 ‰ broad rectangle kernel, and finally interpolated to a 1 ‰ × 1 ‰ grid.Probabilities are normalized to 1 at the maximum; contours show probability levels of 0.8, 0.7, 0.6, 0.45, 0.35.(a) The model assumes sublimation of snow (no isotope fractionation) in the case of moisture uptake and skin temperature (T skin ) below 0 • C. (b) The model assumes evaporation of melted snow (equilibrium isotope fractionation) in the case of moisture uptake and T skin < 0 • C. (c) The model assumes sublimation of snow in the case of moisture uptake and T skin < −7.7 • C. In the case of moisture uptake and −7.7 • C < T skin < 0 • C the model assumes evaporation of melted snow.

Figure 10 .
Figure10.Mean difference between modeled and observed δD of water vapor at Karlsruhe for selected cold snap events ( δD).Each cross represents the δD obtained from one model configuration with specific model assumptions regarding the T subl,max , the δD of snow, the δD at the model initialization, and the amount of moisture uptake.The thick black line connects the δD from 16 model configurations with T subl,max between −15 and 0 • C and the model assumptions of M S_MW,T subl,max and illustrates the increase of the modeled δD values with increasing maximum skin temperature allowing non-fractionating sublimation.The gray shaded area depicts possible systematic changes of the δD in the case of a systematic deviation of the δD of snow from the RCWIP climatology (light blue), a systematically changed δD at the model initialization (yellow), a systematically changed amount of moisture uptake (green), and superposition of the different assumptions (thin black lines).The thin dashed black line corresponds to the model configuration M S_MW,+++,T subl,max , which does not allow reproducing the observations in group "warm" and is therefore not considered in the gray shaded area of uncertainty.In magenta, optimal T subl,max for best agreement of model and observation (dot) and uncertainty of the optimal T subl,max (error bar) due to uncertainty of model assumptions are shown.