Drivers of the variability of the isotopic composition of water vapor in the surface boundary layer

Abstract. Measurements of the isotopic composition of water vapor, δv , as well as measurements of the isotopic composition of evaporation and transpiration provide valuable insights in the hydrological cycle. Here we present measurements of δv in the surface boundary layer (SBL) in combination with eddy covariance (EC) measurements of the isotopic composition of evapotranspiration δET for both δD as well as δO over a full growing season above a managed beech forest in central 5 Germany. Based on direct measurements of isoforcing IF and the height h of the planetary boundary layer (PBL), we provide an estimate of isoforcing-related changes in δv , revealing the influence of local evapotranspiration (ET) on δv . At seasonal time scales we find no evidence for a dominant control of δv by local ET. Rayleigh distillation could at most explain 35 % of the observed variability and we did not find indications for the influence of entrainment at seasonal time scales. Instead, we obtain a strong significant correlation (R ≈ 0.52; p < 10−35) of δv to temperature. We conclude that the observed seasonal 10 variability of δv is neither dominated by Rayleigh processes, entrainment nor local ET but likely linked to other temperaturerelated processes such as fractionation during evaporation. At a diurnal time scale we find that even during summer, when transpiration is high and at a height of only 10 m above the canopy, ET is overruled by entrainment effects throughout the day from approximately 10 am to 4 pm. ET only dominates the diurnal cycle of δv in the mornings and evenings. Thus, from diurnal to seasonal time scale, ET does not dominate the measured δv at our field site, even if the measurements were carried 15 out close to the canopy. We further conclude, that accounting for PBL height h is essential to understand drivers of δv .


Introduction
The isotopic composition of water vapor (δ v ) in the atmosphere can provide insights into the hydrological cycle at scales ranging from the leaf scale to the global scale (Gat, 1996;Jouzel et al., 2000;Yakir and Sternberg, 2000;Wen et al., 2010;Huang and Wen, 2014). Potential drivers of the temporal variability of δ v are the origin and the different histories of air masses 20 (Noone et al., 2013), the removal of rain from the atmosphere, (see e.g. Gat, 2000) as well as horizontal advection and landsea breeze circulation (Lee et al., 2006). Further, δ v is driven by transpiration and evaporation (see e.g. Gat, 2000) and vertical atmospheric mixing, including the growth and decay of the planetary boundary layer (PBL) (Lee et al., 2012a) and in particular entrainment (e.g. Lee et al., 2007). even at a height of 185 m above a cropgrassland canopy, Griffis et al. (2016) estimate the relative contribution of ET to range 55 from 0 to close to 100 %, with a median of 34% (Griffis et al., 2016).
Forest ecosystems have comparably large ET, thus, their contribution to the local variability of δ v is expected to be more pronounced than for other ecosystems. Direct measurements of δ 18 O v and δD v in general and in particular above forest ecosystems are sparse. Only five of the studies in Table 2 (Lee et al., 2007;Welp et al., 2008;Griffis et al., 2010;Huang and Wen, 2014;Griffis et al., 2016) are based on measurements of δ ET and direct eddy covariance (EC) measurements of 60 δ 18 O ET have been carried out only above corn/soybean canopy (Griffis et al., 2010(Griffis et al., , 2011. A forest ecosystem is only captured by one of these studies (Huang and Wen, 2014), based on flux-gradient measurements. Given the limitations of flux gradient  Table 2. Overview of studies that discuss diurnal and seasonal dynamics of δv based on different methods: flux gradient(FG), eddy covariance (EC), repeated profile measurements (RPM), time series analysys (TS) or large eddy simulations (LES). These studies discuss entrainment (ENT), evapotranspiration (ET), dew formation (dew) and land-sea breeze as majos driver of the diurnal variability of δv. The R 2 value of the log linear relationship between δ 18 Ov varies from 0.15 to 0.78. For a list of additional time series measurements of δv, see also (Huang and Wen, 2014

Calculation of isoforcing
We quantify potential effects of ET on the isotopic composition of the atmosphere, by calculating isoforcing, based on EC measurements of the magnitude and the isotopic composition of ET (described in Braden-Behrens et al., 2019). Isoforcing (see e.g. Lee et al., 2009) can be interpreted as the rate of change of the atmospheric δ value multiplied by the boundary layer height h, if a simple isotopic mass balance model (see e.g. Lai et al., 2006) is assumed with only one flux component from the 115 surface and no horizontal advection or entrainment from above (see also Sturm et al., 2012;Braden-Behrens et al., 2019).
With the flux F (e.g. ET), its isotopic composition δ F (e.g. δ ET ), the atmospheric mole fraction C a , the molar density of atmospheric air ρ a , the atmosphere's isotopic composition δ v and the height h of the planetary boundary layer (PBL). shortwave, photosynthetically active radiation and net radiation (see also Braden-Behrens et al., 2019). In this experiment, we also measured the isotopic composition of rain above the canopy using three integrated rain samplers. These rain samplers are self-manufactured according to the description of Gröning et al. (2012). Every second week, we took two samples out of each of the three integrated rain samplers and analyzed them for their δ 18 O and δD composition at the centre for stable isotope research and analysis (KOSI, University of Goettingen) using isotope ratio mass spectrometry (IRMS). After the subsamples were taken, the bottles of the integrated rain samplers were replaced by dry bottles. The uncertainty of our rain sample measurements, as 130 quantified by the median of the standard deviations of the 6 respective subsamples of each data point was approximately 0.2 ‰ for δ 18 O and 1.5 ‰ for δD.

Determination of the PBL height
The time series of the estimate of the height h of the planetary boundary layer (PBL) used in this study is based on a combination of advanced modeling and data assimilation systems delivered by the ERA5 reanalysis product of the European Centre for 135 Medium-Range Weather Forecasts, ECMWF. The height of the PBL was obtained for the full year 2016 at the ERA5 product's native hourly resolution.
The height h of the PBL was extracted from the ERA5 data set at 51.25°latitude, 10.25°longitude, at a distance of 11.94 km from the study site. The grid point used for uncertainty estimation was located at 51.50°latitude, 10.50°longitude, at a distance of 21.18 km from the study site. The spatial grid and the temporal resolution of the reanalysis ensemble uncertainty estimate 140 is coarser (3 h and 0.5 degree resolution) compared to the reanalysis data (1 h and 0.25 degree resolution). Therefore the grid points of the two products do not necessarily coincide, resulting in a different reference grid point for the estimates of PBLheight and uncertainty of PBL-height, respectively, for the current study. In both cases we used the grid point closest to the experimental site.
The following details regarding the PBL data source and processing apply to the current analysis. Here we use the ERA5 145 reanalysis data set generated by 4D-Var data assimilation in CY41R2 of ECMWF's Integrated Forecast System (IFS). The components of IFS are described in detail in ECMWF (2016a,c,d,e,f,g,b only shallow and rather noisy diurnal cycles (Fig. 1). In summer, when the diurnal cycle is most pronounced, δ v is enriched from midnight to approximately 10 am (GMT+1) 5 followed by a depletion throughout the day until 4 pm by approximately 1 ‰ for δ 18 O and 4.5 ‰ for δD.The diurnal cycle of isoforcing values (calculated with Eq. 2) is dominated by the diurnal cycle of ET with a concave shape that is more pronounced in spring and summer and less pronounced in autumn (Fig. 1).
Further, all obtained isoforcing values are positive, thus local ET corresponds to an enrichment of δ v . In summer, the period 165 with the most pronounced diurnal cycles, the amplitude of isoforcing values is approximately 0.08 and 0.6 m‰ s −1 for δ 18 O and δD respectively ( Figure 1). Comparable magnitudes of isoforcing values were found by other authors (Welp et al., 2008;Lee et al., 2007;Hu et al., 2014) in different ecosystems including a temperate forest and a semi-arid area. Diurnal cycles of δ v that are similar to the observed sine-shaped diurnal cycle have been obtained at field sites with less transpiration such as an urban cite in Bejing , or an artificial oasis cropland in the late growing season (Huang and Wen,170 2014). However, a similar diurnal cycle can also be found in a large eddy simulation (Lee et al., 2012a,b). Diurnal cycles of δ v over forest ecosystems have been sparsely measured. However, they have been reported to show constant, concave or convex shapes (Lee et al., 2006(Lee et al., , 2007Welp et al., 2012) and did not always agree between δ 18 O v and δD v (Welp et al., 2012;Lai and Ehleringer, 2011).
The isoforcing-related change in ambient water vapor dδ v /dt iso can be estimated by dividing isoforcing (IF) by the PBL 175 height h (Eq. 2). This estimation is based on assuming a well mixed PBL with constant δ v . 6 Because both, IF and h approach zero at nighttime, the estimated dδ v /dt iso shows large peaks at nighttime. These artifacts are related to dividing by zero. Thus the estimation of dδ v /dt iso only provides meaningful values at daytime. Throughout the day, from 8 am to 5 pm, the estimated dδ v /dt iso fluctuates around 0.1 permil for δ 18 O and 1.1 permil for δD, with slightly lower values in autumn (see Fig. 2).
The directly measured diurnal cycles of dδ v /dt meas do not agree with the isoforcing-related estimate dδ v /dt iso (see Fig. 2). In 180 particular in spring and summer, we measure negative values of dδ v /dt around midday, associated with a depletion of ambient water vapor, while the isoforcing-related change δ v always yields an enrichment. However, the magnitude of the diurnal cycle of dδ v /dt meas is comparable to the isoforcing-related estimate dδ v /dt iso .
The mean diurnal cycles of our dataset ( Fig. 1 and Fig. 2) show conclusively that isoforcing due to local ET does not dominate the diurnal variability of δD v and δ 18 O v . This is also the case during summer, when transpiration is high -even if our 185 4 https://confluence.ecmwf.int/pages/viewpage.action?pageId=111158117 5 All times given in this paper refer to local winter time (GMT+1). 6 This assumption is probably not justified, but the calculation of dδv/dt should provide a rough estimate of the influence of local ET on δv.   measurements took place only 10 m above the tree tops. Thus, the measured data do not support the hypothesis that during summer, local ET dominates the diurnal variability of δ v over a forest ecosystem. We further conclude that due to the large variability of the boundary layer height h, it is essential to account for h when estimating the influence of local ET on ambient water vapor. A discussion of the influence of local ET that is purely based on isoforcing IF overlooks the influence of boundary layer mixing processes. Further, the concurrent trends in the diurnal cycles of C H2O and δ v (Fig. 1) indicate, that entrainment 190 dominantly influences δ v from the forenoon to the afternoon: The rise in δ v in the morning and in the afternoon is related to a water source (e.g. local ET) whereas the decrease of δ v throughout the day is related to mixing with dryer and isotopically lighter air (entrainment). Despite our measurement position beeing closely above the forest canopy and despite high transpiration, we observe this indication for a dominant influence of entrainment from the forenoon to the afternoon also in summer.
nately dominated by entrainment and local ET. As turbulence increases boundary layer mixing, we suppose that larger TKE corresponds to larger entrainment. This is also supported by the strong agreement between the diurnal cycles of TKE and PBL height h (Fig. 2). They co-vary throughout the day with an increase in the mornings, when the PBL grows, and a maximum around 2 pm. The observed negative values for dδ v /dt iso throughout the day occur during a time period with large TKE (Fig.   2). Additionally, the amplitude of TKE is largest in summer, medium large in spring and comparably small in autumn (Fig. 2). δ v and the corresponding isoforcing values IF as well as the isoforcing related change dδ v /dt iso (R 2 ≈ 0.2) for both, δD v and δ 18 O v . We propose that these weak correlations between isoforcing-related quantities do not reflect a causal relationship.
Additionally to being weak, they are all negative (as shown in Fig. 5). If local ET drove the isotopic composition of ambient water vapor, we would obtain a positive correlation between IF and δ v . Thus, we discard the hypothesis that at the field site of our measurement campaign ET dominates the seasonal variability of δ v . We further conclude, that it is important to use direct 220 measurements of IF and an estimation of dδ v /dt iso to not over-interpret correlations between δ v and δ ET . In our dataset, we indeed find a positive moderate (R 2 ≈ 0.4) and significant (p < 10 −25 ) correlation between δD ET and δD v , but not for δ 18 O ET and δ 18 O v (see Table 3).
Rayleigh rain-out is another potentially important driver of the seasonal variability of δ v . Rayleigh distillation processes yield a log-linear relationship between δ v and C H2O (Eq. 1) as tested in Fig. 6. Throughout the whole measurement period, 225 both δ v values correlate only moderately (R 2 ≈ 0.35), but significantly (p < 10 −25 ) to log(C H2O ). This correlation is weaker 7 Here we choose the midday value of dδv/dt iso dt to reflect dδv/dt iso dt. The reason for this choice is that due to data gaps and due to the peaks in dδv/dt iso at nighttime, the numeric integration of dδv/dt iso would not be meaningful. in summer (see Table 3). Stronger correlations between δ v and log(C H2O ) at the seasonal time scale have been reported for various other sites, for δ 18 O (Lee et al., 2006(Lee et al., , 2007Welp et al., 2008;Wen et al., 2010;Zhang et al., 2011;Griffis et al., 2016) and also for δD (Wen et al., 2010;Zhang et al., 2011). Similarly small correlation coefficients (R 2 < 0.17) that get weaker during summer, were obtained above an arid artificial oasis (Huang and Wen, 2014). At our field site, the correlation between 230 δ v and log(C H2O ) is dominated by the periods before leaf unfolding and after leaf senescence (Fig. 6). During the period between leaf unfolding and senescence, this correlation is particularly weak, with R 2 values below 0.2 (Table 3). Thus, our  before leaf unfolding (orange crosses), the period with green unfolded leaves (green triangles) and the period after the beginning of leaf senescence in fall (brown diagonal crosses). A log-linear relationship would indicate a system that is dominated by Rayleigh distillation.
We plot the measured δ v values in the δ 18 O-δD-plane (Fig. 7), to evaluate deviations from Rayleigh distillation processes.

235
Therefore, we analyze deviations of measured δ rain and δ v from the global meteoric water line 8 (GMWL, δD=8δ 18 O+10) (Craig, 1961;Dansgaard, 1964;Gat, 2000;Bowling et al., 2017), that corresponds to equilibrium fractionation during Rayleigh rain-out (Dansgaard, 1964). At our field site, the the local meteoric water line (LMWL), as defined by the linear regression of rain sample data (Dansgaard, 1964;Bowling et al., 2017), has a slope of approximately 7.4 ± 0.3. This lower slope of the LMWL compared to the GMWL might reflect the influence of non-Rayleigh-distillation processes such as local evaporation

240
(from open water bodies) and selective transpiration but also evaporation from falling raindrops (Gat, 2000). However, the LMWL is to be interpreted cautiously, as the underlying rain samples span different seasons (see Gat, 1996). Here we use the 8 The GMWL is based on meteoric waters, i.e. precipitation and waters from rivers and lakes (except East African rivers and lakes) and waters from closed basins (Craig, 1961). LMWL and the GMWL only to compare the slopes of the measured isotopic composition of water vapor with these lines (Fig.   7). The dual isotope analysis reveals that the measured δ v values over the season formed the clusters: Before leaf unfolding on DOY 110, the measured (daily averaged) δ v values followed the GMWL (with a slope of 8.0 ± 0.2), whereas after leaf 245 senescence on DOY 280, the slope in the δ 18 O-δD-plane was 7.3 ± 0.1 and thus closer to the LMWL. For the period with green leaves, an even lower slope of 6.9 ± 0.1 was measured. This becomes more visible when deviations from the GMWL and LMWL are plotted versus time (Fig. 7, right panels). While the GMWL with its slope of 8 seems to describe the variability of δ v before leave unfolding, the obtained values after leaf senescence in fall are on average better represented by the LMWL, indicating some influence of local conditions or fluxes. Between leaf unfolding and leaf senescence in fall, the variation in the 250 measured difference from both the LMWL and the GMWL shows a seasonal cycle which might be related to seasonal shifts in the source (but in general also in the fractionation) of water vapor.
To reveal other processes that drive the seasonal variability of δ v , we calculate the Pearson correlation coefficient R pear for different meteorological quantities (Table 3). We find no indication that the seasonal variability of δ v is driven by the variability of entrainment, as we find no significant correlation between δ v and turbulent kinetic energy TKE and friction velocity u * (see 255   Table 3). For the whole measurement period, the observed seasonal variability of δ v was strongly correlated to temperature (R 2 > 0.5 , p < 10 −35 ). This correlation was stronger than the correlation to log(C H2O ), discussed above as an indicator for Rayleigh distillation processes. In particular, the observed correlation with temperature at 2 m height above the ground had an R 2 of 0.52 and was slightly stronger than the correlation with temperature at 44 m above the ground for both isotopic species.
Similarly, soil temperature at 2 cm depth is much stronger correlated to δ v than soil temperature at 64 cm depth. These height 260 Table 3. Results of the analysis of potential drivers of δv for all data points ('all times') and for the period with green unfolded leaves ('green leaves'). For each correlation the sign and the R 2 value is given in combination with the significance levels, marked with (*) for p < 10 −5 , (*) for p < 10 −10 , (*) for p < 10 −15 and so on. dependencies indicate, that the temperature close to the surface is an important driver of δ v . In general, the correlation between temperature and δ v might be linked to temperature dependent fractionation at the sites of evaporation. However, the day-to-dayvariability is not fully reflected by the obtained correlations of δ v with temperature and temperature-related quantities such as period, the obtained correlations with temperature-related quantities get weaker. The correlation with temperature at 2 m height for example is still significant (p < 10 −10 ) but has R 2 -values of only 0.27 and 0.28 for δD v and δ 18 O v , respectively. In general, the positive correlation with temperature related quantities implies that fractionation and evaporation might be relevant drivers of δ v . As we did not find indications that local ET drives the seasonal variability of δ v , the dominant source of the measured water vapor might be further away, but the temperature during evaporation might still be correlated to the temperature during 270 the measurement.

Conclusions
Here we evaluate laser spectroscopic measurements of the isotopic composition of water vapor (δD v and δ 18 O v ) in the SBL and local ET (δD ET and δ 18 O ET ) in combination with meteorologic and turbulence-related quantities on diurnal to seasonal time scale above a managed beech forest in central Germany. We find that the temporal variability of the isotopic composition 275 of water vapor in the SBL, even if measured closely above a managed beech forest, is not dominated by local ET at both, diurnal and seasonal time scale. Based on directly measured δ ET and PBL height h extracted from reanalysis data, we present an estimate of isoforcing-related change in the isotopic composition of ambient water vapor, dδ v /dt iso . A direct comparison of measured dδ v /dt with isoforcing-related dδ v /dt reveals that local ET does not dominate the diurnal cycle of δ v throughout the year. We conclude that the diurnal cycle of δ v is alternately dominated by local ET and entrainment, in particular in 280 spring and summer. At seasonal time scales, we find no indication that local ET and entrainment dominate the observed variability. This could be related to a mutual canceling out of entrainment and local ET at this field site, yielding a more complicated seasonal variability. Further, the variability of δ v over the full growing season can only partly be explained by Rayleigh distillation (linked to approximately 35 % of the variability). This fraction further decreases for the period when green leaves were present and transpiration is expected. A larger fraction of 50 % of the observed seasonal variability of δ v 285 is linked to temperature, indicating some influence of temperature-related processes, such as fractionation during evaporation from different water pools. We conclude that the simultaneous measurement of δ v and δ ET in combination with meteorological and turbulence-related quantities and PBL height h seems a promising approach to increase the understanding of the temporal variability in δ v . As the influence of local ET on δ v is strongly moderated by boundary layer processes, accounting for PBL height is particularly essential to understand changes in δ v . uncertainty of less than 31% of the PBL height h, with lower relative uncertainties for higher PBL heights (see Fig A1) in the appendix. The mean uncertainty of the PBL height h during 2016 was 6.6% of PBL height h, obtained as the slope of a major axis linear regression of the uncertainty of h versus h itself. Thus, the average uncertainty of 6.6% is small compared to the typically large diurnal changes of the boundary layer height. Given above uncertainty estimates of the reanalysis data, it is appropriate to estimate boundary layer height based on reanalysis data.