Determining water sources in the boundary layer from tall tower profiles of water vapor and surface water isotope ratios after a snowstorm in Colorado

Introduction Conclusions References


Introduction
The representation of the land surface energy and water budget is a significant source of dispersion between climate model projections of future changes in temperature and precipitation in response to greenhouse forcing (e.g., Crossley et al., 2000;Gedney et al., 2000;Boe and Terray, 2008) or land use changes (Pitman et al., 2009).Land-atmosphere interactions also play a key role in the meteorological and hydrological variability at different time scales from weeks to decades (Avissar and Werth, 2005;Milly et al., 2005;Bosilovich and Chern, 2006;Guo et al., 2006;Koster et al., 2006).Inadequate representation of surface fluxes and their dependence on surface conditions are among the key sources of uncertainties in quantifying regional hydroclimate (Koster and Milly, 1997;Gedney et al., 2000;Boone et al., 2004).

D. Noone et al.: Water sources in the boundary layer
A particularly vexing subset of this issue arises in regions covered by snow, where the representation of processes in the snow pack is largely unconstrained (Slater et al., 2001;Boone et al., 2004;Rutter et al., 2009).Additional surface flux uncertainties arise from imperfect depiction of boundary layer processes, notably under stable conditions (Holtslag, 2006).Reducing these uncertainties requires better understanding of the movement of water between the landscape and the atmosphere.
Because the stable hydrogen and oxygen isotope ratio of water vapor, liquid and ice reflects the balance of processes influencing regional hydrology, measurements of the oxygen and hydrogen isotope ratios of water can provide constraints for water balance studies and expose model shortcomings (Jouzel and Merlivat, 1984;Henderson-Sellers et al., 2004;Henderson-Sellers, 2006;Sturm et al., 2010).The utility of isotope ratio information stems from the fractionation that accompanies phase changes in which the heavier isotopologues preferentially remain in liquid or solid form during evaporation and condensation (Bigeleisen, 1961;Dansgaard, 1964).At the continental scale, the distributions of the isotopic composition of precipitation have been used to partition continental recycling into evaporation from standing water and transpiration (Salati et al., 1979;Gat and Matsui, 1991).At the local scale, isotopic measurements of water vapor and soil water have been used to partition evapotranspiration into transpiration from plants and evaporation from soil (Moreira et al., 1997;Yepez et al., 2003;Williams et al., 2004;Wang et al., 2010).Several water vapor isotopic measurements at the surface or within the boundary layer have suggested that the isotopic composition could also provide information on turbulent properties (Ehhalt, 1974;He and Smith, 1999;Ehhalt et al., 2005;Angert et al., 2008;Noone et al., 2011).However, studies to date have been limited by the availability of measurements of water vapor isotopic ratios.
A common approach in water source partitioning is to estimate the isotopic composition of the total surface flux.This is typically estimated using a " Keeling plot" approach (Keeling, 1958;Pataki et al., 2003) which is based on a simple two-member mixing model.This technique has been applied to both temporal series (e.g., Moreira et al., 1997;Noone et al., 2011;Noone, 2012) and vertical profiles (Yepez et al., 2003;Williams et al., 2004).However there is some evidence from CO 2 and 13 C measurements from tall towers that both have limitations due to advection that invalidates the assumptions underlying the method (Lee et al., 2006(Lee et al., , 2012;;Griffis et al., 2007).With the advent of field deployable water vapor isotopic instruments, many of the previous technical limitations in obtaining water vapor isotopic data have been overcome and high frequency measurements can now be performed (e.g., Lee et al., 2007;Gupta et al., 2009;Wang et al., 2010).Such measurements allow reexamination of the mixing line method and the question of how to attribute fluxes to multiple contributing components.It is anticipated that gradient theory over the depth of the boundary layer will not provide unambiguous estimates of surface fluxes because of the non-stationary nature of the boundary layer flow.There remains a need to identify the depth over which estimates of the surface layer fluxes are valid versus the need to account for the more dynamic flow regime occurring higher.
In this paper, we evaluate some features of the boundary layer relevant for understanding the transport of water and trace gases.We test the use of isotopic information to decompose surface latent heat flux into different physical components.The flux attribution calculation employs mass balance to partition the flux associated with loss of the snow pack into sublimation and evaporation from melt ponds and employs a maximum likelihood estimation approach to provide an optimal estimate in the under-constrained problem.The technological improvement in observational capacity is exploited to measure vertical profiles of water vapor isotopic composition within the boundary layer from the surface to 300 m near Boulder, Colorado at a vertical resolution of a few tens of meters every 15 min.The study focuses on results from a four-day field campaign in February 2010 during which a fresh surface snow pack of about 10 cm depth vanished.Before undertaking the mass balance estimate, the data are used to evaluate the applicability of the mixing line method in investigating the controls on water vapor isotopic composition on diurnal and synoptic time scales.We conclude that with greater amounts of water vapor isotopic information, limitations of mixing line methods are exposed and more advanced models and estimation techniques are warranted.

Measurements and data
A four day field campaign was held between 15 and 18 February 2010 at the 300 m turbulence research tower at the NOAA Boulder Atmospheric Observatory (40.050 • N, 105.003 • W, 1584 m a.s.l.).The location is mostly flat terrain about 25 km east of the foothills of the Rocky Mountains (Kaimal and Gaynor, 1983).The facility has been used for a large range of atmospheric research applications such as boundary layer meteorology (Blumen, 1984;Gossard et al., 1985), wave-turbulence interactions (Einaudi et al., 1989;Einaudi and Finnigan, 1993), atmospheric chemistry (Brown et al., 2007), and instrument testing (Cohn et al., 2001).The footprint for flux measurements is typically around 30-50 km, estimated from both turbulence theory and from known industrial and vehicle emissions from nearby cities, including Boulder (P.Blanken, personnel communication, 2012).The February period was chosen because it followed a winter snow storm that ended on 14 February, and the experiment terminated after nearly all the snow had been lost from the surface and snow had begun to fall in association with another storm on the evening of 18 February.The weak synoptic evolution and moderate local wind speeds between the snowstorm of the 13-14 February and the subsequent snow storm on the 18 February provided an ideal opportunity for the boundary layer and surface water balance study.
During the experiment, a water vapor isotopic analyzer (Picarro model L1115-i, Gupta et al., 2009) was installed on an instrument elevator platform and measurements were made continuously at approximately 0.16 Hz.The elevator platform moves the length of the tower over a nineminute period with an ascent/descent rate of approximately 0.55 m s −1 .During the experiment, profile measurements were started every 15 min (nine minutes for the elevator to ascend or descend on the tower, and six minutes when the elevator was stationary) to give a total of 312 profiles in either the "up" or "down" direction.Time series measurements were converted to height profiles using heights derived from the (assumed constant) ascent/descent rate.The isotopic analyzer was installed in a non-temperature controlled enclosure with a heated 0.125 inch outer diameter stainless steel sample inlet line drawing air into the analyzer at 30 scc min −1 mass flow from a horizontal boom that extended approximately 4 m from the tower.The inlet tip was the top half of a plastic bottle which was used to prevent any precipitation from entering the sample line.Isotope ratios, R, are reported in "delta" notation (δ = R/R std − 1, and R std is the Vienna Standard Mean Ocean Water standard), and δ 18 O is for R = 18 O/ 16 O and δD is for R = D/H.The raw isotopic δ values were corrected to sequentially account for (1) humidity dependence, (2) calibration to the primary reference scale and (3) memory effects.Careful calibration is needed for spectroscopic measurements, and some details of the approach used here are given in Appendix A. Memory in the measurement system can arise from the instrument's internal plumbing, in association with wake turbulence near the boom arm and the gas volume in the inlet and sample line.The algorithm to minimize the influence of memory effects is given in Appendix B. There were no gaps in the data over the duration of the experiment, and there appeared no need to reject data based on wind direction due to interference from the tower.The average accuracy of corrected values used in the analysis is 3.4 ‰ for δD and 0.45 ‰ for δ 18 O.Other groups report similar performance using similar instruments (e.g., Gupta et al., 2009;Noone et al., 2011;Rambo et al., 2011;Kurita et al., 2012;Tremoy et al., 2012;Wen et al., 2012).There is significant potential for sensitivity for temperature variations during diurnal cycles and instrument vibration to limit the instrument performance.This has been demonstrated, for instance, in some field settings (Aemisegger et al., 2012).Rather than being resolved and corrected, drift effects becomes part of the total uncertainty estimate.We find that these are particularly important for deuterium excess, rather than either single isotope ratio δD or δ 18 O.Due to this potential limitation we focus here only on a single isotope ratio, δD, and use the deuterium excess information only as a check on the δD value and to ensure physical consistency in conclusions.
Surface snow, puddle, mud and frost samples were collected into 60 mL wide mouth Nalgene bottles and stored frozen until time of analysis in our laboratory.Water samples were obtained from mud using the cryogenic vacuum extraction method based on published protocols (West et al., 2006), with two repeat extractions from each mud sample yielding around 1 mL of water.Samples were weighed before and after extractions to ensure the water mass transfer was complete.All liquid isotopic analyses were performed using the same L1115-i isotopic analyzer that was used on the tower.The instrument precision for liquid analyses is less than 0.5 ‰ for δD and 0.02 ‰ for δ 18 O.Each sample was injected three times to allow correction for memory effects.Instrument drift was corrected using injections of a standard after every eight analyses, and calibrated to the absolute scale using a three-point calibration (Florida Water, Boulder Water and West Antarctic Divide secondary standards tied to the International Atomic Energy Agency scale with mass spectrometer determinations, B. Vaughn and J. White, personal communication, 2010).Although measurements of both δ 18 O and δD were made, the analysis here focuses on the single species δD.
Profile measurements of temperature and wind speed were made using a very fast response (response time < 0.001 s, and sampled at 1000 Hz) thermocouple and hot wire probe on the Cooperative Institute for Research in Environmental Sciences tethered lift system (TLS) to characterize the turbulence structure of the lower boundary layer (Balsley et al., 1998;Frehlich et al., 2003).The TLS comprises the instrument package suspended from a helium filled blimp with profiles of temperature and wind speed attained during slow ascent and decent (0.1-1 m s −1 ) that was digitally controlled by a winch.These measurements were made approximately 400 m to the southwest of the main tower, and only data from the night of 16 February are reported here.
Additional measurements from the tower elevator included high speed (10 Hz) wind from a Gill Windmaster ultrasonic anemometer and H 2 O and CO 2 concentrations from an open path non-dispersive infrared analyzer (Licor 7500) provided by the NOAA Physical Sciences Division.The H 2 O mixing ratio measurements were in excellent agreement with the Picarro measurements once calibrated, as anticipated from the previous work (Noone et al., 2011).Temperature measurements on the platform are from a Vaisala HMP 45C sensor, which has a characteristic response time of less than one minute.While vibrations of the elevator are likely to influence the high frequency (faster than 1 Hz) measurements, our analysis focuses on 10 s averages to minimize vibration influences.There were no surface sensible or latent heat flux measurements made during the period of the study, while radiation data show characteristically low energy availability due to the high albedo of snow.Noone et al. (2011) showed that representing turbulent exchange as a diffusive process leads to mixing models that can be applied to isotope ratios as either spatial gradients or time series.The result can be derived simply from the mass balance for a single atmospheric parcel.Specifically, the water vapor mixing ratio (q) is considered the sum of some background water vapor (say, free troposphere water vapor mixing ratio, q T ) and water vapor derived from a surface flux (q F , either a source or sink) such that: (1) A similar mass balance can be written for mixing ratio of HDO.Using equation 1 to eliminate q F and converting to δD one can write where again the subscripts T and F denote the free troposphere and the surface flux.This is analogous to the wellknown derivation given by Keeling (1958), and suggests a graphical method determining δD F by plotting measured δD as a function of the reciprocal of q and finding the intercept (i.e., the asymptotic limit at q −1 = 0).While the q −1 plot is a convenient graphical device, Miller and Tans (2003) showed that regression errors are reduced if Eq. ( 2) is written as and δD F is found as the slope of the regression of qδD versus q, with error in δD F described by the confidence limits on the regression coefficient given total uncertainties (quadrature sum of accuracy and precision) in both q and δD.This model can be applied freely when the vapor flux and the free tropospheric vapor remains unchanged with time.The model may also be valid in some instances even when the tropospheric background is changing (Miller and Tans, 2003).The surface flux can be from evapotranspiration or sublimation and Eq. ( 2) is also applicable to frost (or dew) with the sign convention that q F is negative.In this case, δD F gives the isotopic composition of the frost.Mixing line analysis is often applied to time series data [hereafter identified as the temporal mixing line (TML) method], and requires that time variations in δD are due to the progressive input of the same water vapor source.This is not always a good assumption.
It is useful to note that a mixing line analysis can be reconciled with similarity theory within the surface layer.The steady state profile of wind and constituents (such as H 2 O, HDO, CO 2 ) in the constant flux surface layer is and where z is height, z 0 is the roughness length, k is the von Karman constant (∼ 0.4), q 0 is the mixing ratio analogous to z 0 , L is the Obukhov length, is the turbulent structure function for momentum (subscript m) or trace constituents (such as water, subscript q) from similarity under non-neutral conditions, u * is the friction velocity and q * is the moisture perturbation scale.Given an air density of ρ, the latent heat flux (either positive or negative) is E = −ρu * q * .For the isotopologue mixing ratio the comparable relationship is where q may be used for all isotopologues when turbulence dominates.Analogously, the evaporative flux of the isotopologue is E i = −ρu * R * q * , and R * is the isotope ratio of the flux.As with u * and q * , R * can be readily evaluated as the slope of a linear regression of observed values of Rq as a function of the term in square brackets [≈ ln(z)].We refer to this as the gradient mixing-line (GML) approach.
Equivalence between the GML and the one-box temporal mixing model (i.e., TML given by Eq. 2) is immediate by substituting equation 5 into 6 to eliminate the height term and q * , dividing by q then subtracting 1 to convert the isotope ratios to δ values to obtain: This result suggests an isotopic mixing line approach is appropriate in the surface layer when the mixing length hypothesis is valid.It suggests the profile measurements (e.g., Yepez et al., 2003;Williams et al., 2004) may be more generally applicable than the TML method which is based on time evolution of a single control volume as in Eq. ( 2) (e.g., Moreira et al., 1997) and requires more ideal conditions to be valid (e.g., Good et al., 2012).The validity of the time versus space depiction of mixing requires testing with field data.Note that similar relationships between air mass mixing derived from spatial or temporal data, including the special case of the surface layer, are explained by the general solution examined elsewhere (i.e., Noone et al., 2011).
Given that the simple mixing configurations assume stationary end members, assessing the fidelity of mixing line methods comprises identifying changes in the mixing line end members.For profile data obtained we consider three distinct configurations where a single factor dominates the measurements at any time (Fig. 1).First, a surface source associated with evaporation or sublimation will increase the humidity and, likely, the isotope ratio from the bottom up (Fig. 1a).Second, in the case of near-surface formation of dew or frost, the mixing is between the ambient vapor and the vapor that remains after the continual loss of water and heavy isotopes from the surface vapor (Fig. 1b).Third, in the absence of surface sources, vertical mixing at the top of the boundary layer with the free troposphere at night leads to mixing between the initial boundary layer vapor and the free troposphere vapor (Fig. 1c).These situations highlight that the mixing approach is built on the assumption that the profile structure is not influenced by other factors such as lateral inhomogeneity and advection, which we show confounds the technique.Indeed, previous work on boundary layer flow and trace gas transport has illustrated it is unlikely that the simple mixing model is valid for many instances in the lower boundary layer to, for instance, entrainment and lateral advection.

Mass balance for surface water and fluxes
Estimates of δD F are used to attribute the surface water flux into contributions from different components.A mass balance model for surface snow and liquid is depicted in Fig. 2 and is motivated by visual evidence for changes in snow grain size and the formation of muddy ponds during the experiment.The quantitative mass balance assumes that over some time interval a portion of the snow pack melts (f melt ) and a portion sublimates (f sub ).Some portion of the melt water in the snow drains and adds to the mass of ponds (f drain ) and a portion can evaporate (f evapsnow ).Similarly, some portion of the pond water evaporates (f evappond ).The model does not account for water drainage into the soil, which is likely small given the dry frozen clay soils at the site.The fractional values of these portions are constrained by observations of the snow and pond water isotope ratios, the estimate of δD F obtained by a mixing line analysis described in the previous section, and a model to describe the isotopic exchange be- Schematic depiction of the surface water balance in a region of snow pack loss.Total evaporation is the sum of sublimation, evaporation from melt water in the snow, and evaporation from melt water in ponds on the soil.It is assumed that additional evaporation from preexisting soil water can be neglected.Melted snow may remain in the snow pack or drain to (muddy) ponds.Isotopic fractionation occurs during exchange between the snow pack liquid and ice, and during evaporation from liquid pools.
tween the snow, melt water, pond water and the ambient vapor.
An enrichment of snow (by about 20 ‰) and ponds (by about 35 ‰) was observed at the end of each day compared to morning snow.The enrichment of pond water is likely due to evaporation.The enrichment of snow water can have several causes.First, isotopic re-equilibration between snow crystals and melt water as it percolates through the snow pack can lead to a significant enrichment (Taylor et al., 2001;Lee et al., 2010).The isotope ratio of ice, R c , in the case of partial re-equilibration is where R l0 and R c0 are the initial isotope ratios of the liquid and ice, α f is the liquid/ice fractionation coefficient (ice is 19.5 ‰ enriched compared to liquid at equilibrium, O'Neil, 1968), a is the mass of ice affected by the equilibration per unit mass of liquid.Knowing the time evolution of R c , a can be deduced.Second, snow enrichment can occur due to recrystallization of melt water within the snow pack that has undergone partial evaporation (Gurney and Lawrence, 2004).Third, there remains the possibility of fractionation during snow sublimation, (Ekaykin et al., 2009), which may incur a kinetic effect and would lead to enrichment as in a Rayleigh process.We neglect this possible effect in the default configuration of the mass balance model because the process is not well constrained by existing theory and because it can be shown that it is not necessary to explain the enrichment.increase heterogeneity, but is likely to be a small contribution to the mass.By design, snow samples were collected from a single depth during the study to provide a bulk characterization of the isotope ratios in the diminishing snow pack.While heterogeneity and layering of isotope ratios within the snow and on the landscape is likely important, this aspect cannot be checked with the samples available.
The isotope ratio of evaporating standing water can be described following Stewart (1975).Assuming the reservoir of boundary layer vapor (R v ) exchanging with ponds is much larger than the reservoir of pond water, the isotope ratio of the pond liquid (R l ) is governed by the combination of a Rayleigh-like distillation and re-equilibration of the liquid with the ambient vapor.The model is where α l is the liquid/vapor equilibrium fractionation coefficient (Majoube, 1971), α k is the kinetic fractionation coefficient, h is the relative humidity, 1−f evappond is the fraction of pond water that remains after evaporation, and R l0 is the initial liquid isotope ratio.The kinetic fraction is calculated as where D and D i are the diffusivities of light and heavier isotopic species in air respectively (Merlivat, 1978) and n is a coefficient ranging from 0.58 (in the case of standing water, Stewart, 1975) to 0.67 (in the case of saturated soil, Mathieu and Bariac, 1996).The value of n parameterizes the influence of the balance between turbulent and molecular transport on the isotope ratios.While this treatment of kinetic effects is common in the literature, it is incomplete since it does not include variations in the laminar boundary layers and roughness at the ground or at the surface of snow grains.Nonetheless, here we use n = 0.58, and testing shows this choice has little effect on the results for a single isotope (i.e., δD).The influence of this choice on the deuterium excess is larger, but not discussed here.The model has six unknown parameters: f sub , f evapmelt , f evappond , f melt and f drain , and the parameter a in Eq. ( 10).They are obtained by an optimal estimation approach that is constrained with observations of R v , R l , R c , and the estimate of R F derived from the mixing line analysis.Since the snow pack and nearly all ponds disappeared over 3 days we assume a third of the mass was lost on each day between 06:00 a.m. and 06:00 p.m.The rate of loss of snow mass is assumed constant with time.It is assumed that at 06:00 a.m., all the snow mass loss is due to sublimation given that cold morning temperatures preclude standing liquid, but the proportion of snow sublimation versus melt decreases linearly to a fraction that is two times f sub (the factor of two is so the daily mean sublimated fraction is f sub ).The model is integrated 10 000 times using uniformly randomly chosen parameters that sample their full possible range.Simulations that yield predictions within the (one standard deviation) uncertainty of the observations are selected as plausible solutions.The mean and spread of the ensemble of plausible outcomes are reported as the optimal set which describes the flux partitioning.This optimal estimation approach represents a significant conceptual advance in methods that seek to use isotopic information as a constraint on mass budgets.

Overview of synoptic evolution
A snow storm on 12 February (36 h before the campaign) was accompanied by a cold air mass and synoptic-scale flow from the northwest.Flow over Colorado remained mostly north-easterly at upper-levels between 14-17 February with weak anticyclonic flow in the wake of the earlier front.Subfreezing temperatures of −12 • C during the night on 15 February were associated with both upper level cold air advection from the north (not shown) and total near-surface air temperature decreased by 15 • C between local sunset and the early morning temperature minimum.During the day on 18 February, upper level clouds were associated with westerly flow while high humidity air at lower altitudes was associated with southerly advection.Fresh snowfall began at around 05:00 p.m. local time on 18 February and ended the experiment.Figure 3 shows time-height cross sections of temperature, wind speed, water mixing ratio (q), δD and CO 2 mixing ratio.The water mixing ratio is from the Picarro analyser for consistency with the δD data, however it agrees well with the higher frequency open path Licor measurements.In part the agreement comes from robust calibration of the two instruments; however, the agreement also provides confidence in the memory correction used during the calibration of the Picarro water and isotope ratio data.Time series of specific humidity, δD and CO 2 mixing ratio at the surface and at 300 m are shown in Fig. 4.
At 300 m, both specific humidity and δD increase during the afternoon on 17 February (Fig. 4), and signifies a change in the air mass.The isotopic signature of this transition is characterized using the position of the observed q, δD pair on a q − δ diagram (Fig. 5).This type of diagram helps identify the processes controlling the isotopic ratios (Noone, 2012).The data provide evidence of two distinct moisture pathways.The difference in the two populations of observations is expected from the change in the synoptic flow to the west and south on 18 February.For most of the experiment, the values remain close to a single coherent area of the phase diagram clustered around a line that is consistent with a saturated Rayleigh distillation (i.e., a pseudoadiabatic process) from an oceanic moisture source at 10 • C and 70 % relative humidity (purple curve).At the end of the experiment the air mass associated with the snowfall on the 18 February lies on a Rayleigh curve drawn with a saturated ocean source at around 25 • C (solid cyan curve).A model in which the precipitation efficiency is 0.5 (i.e., partial removal of condensate shown as dashed cyan curve, see Noone, 2012) captures the values early in the onset of the storm, suggesting a history in which liquid cloud processes or recycling are important during transit from the ultimate oceanic source.The distillation curves on the figure were visually fit to the data to provide guidance on understanding origins of the air from the raw observations.Simple single-process modeling approaches can be confounded by more complex air mass histories involving, for instance, air mass mixing.Air mass mixing during transport can lead to substantial deviations from a Rayleigh prediction (e.g., Noone, 2008;Hurley et al., 2012).Indeed, the transition between the two air masses is well described by a mixing line constrained on the upper end with an end A further change in air mass occurred during the snow storm on 18 February that ended the experiment.The shift has two components evident in the δD observations at the end of the time series.First, enrichment is associated with the air mass change to the westerly flow in the morning (as noted above).Then, a strong depletion is associated with strong vertical mixing involving the entire troposphere in frontal systems and is accompanied by high wind speed in the final hours of the experiment (Fig. 3).Values observed during the event are well described by the Rayleigh curve (solid cyan, Fig. 5), and indicates a continually precipitating air mass with an oceanic origin in the region of the front.

Diurnal evolution
Many salient features observed are related to the diurnal cycle.Figure 4 shows humidity at both the surface and 300 m decrease linearly throughout the first two nights, accompanied by a similar steady linear decrease of δD values.On the first night this is associated with stable subsidence and the second night this is associated with the formation of a low level jet (LLJ).At night, surface CO 2 mixing ratio increase with time due to continual respiration while 300 m values remain rather constant and are indicative of stable stratification.CO 2 builds slowly in the stable nocturnal surface layer to a height of about 50 m, and is shown to be approximately bounded by contours of potential temperature in Fig. 3.In the morning, specific humidity and δD values increase due to the surface latent heat flux and the growth of a convective boundary layer which allows ventilation of near-surface air with low CO 2 and low δD air entrained into the boundary layer from the troposphere.The convective boundary layer is Theoretical curves are placed subjectively to provide context.The purple curve is a saturated Rayleigh curve assuming source vapor is initially at 70 % relative humidity over an ocean at 10 • C as might be typical of a North Pacific origin, and the cyan line shows a saturated Rayleigh curve assuming the source vapor is saturated over an ocean at 25 • C as is more typical of a Gulf of Mexico origin.The dashed cyan curve has the same origin, but the precipitation efficiency is 0.5.The orange line is for mixing between an air mass with q = 0.8 g kg −1 and δD = −350 ‰ and a source with δD = −190.
shown on the figure by the height at which the bulk Richardson number attains the critical value, and indicates deepening during each day due to shear.The convective layer exceeds 300 m on the last day.In the evening, water vapor mixing ratio and δD decreases at 300 m once thermally driven convection terminates which stops the steady input of surface water of high δD.The surface source of high δD vapor does not cease immediately in the evening but becomes trapped within the growing stable boundary layer and does not reach 300 m.
Two-hour averaged δD profiles are shown in Fig. 6 for different times of day that illustrate the evolution of the profile in the cases anticipated from Fig. 1.In the morning, sublimation of snow increases the humidity and δD in the boundary layer.During the night the boundary layer attains lower δD values associated with air entering from above (Fig. 6a).Evidence for an evaporative source is seen in the δD as a maximum near the surface (Fig. 6b).This moist and enriched anomaly propagates quickly upward due to the strong mixing that accompanies the steadily growing convective boundary layer and leading to a shift of the δD profile to higher values (Figs. 3 and 6b).The CO 2 , water vapor and δD all show reduced vertical gradients during the day associated with strong mixing driven by solar heating (Fig. 3).In the meantime, mixing within the boundary layer and with the free troposphere allows the night time CO 2 maximum to dissipate.

Evolution of turbulence properties
Since high δD water vapor is emitted at the surface during the day and CO 2 is emitted at the surface is trapped during the night, δD and CO 2 are complementary tracers that reveal information about boundary layer mixing processes.For stationary boundary layers, the relative role of mixing within the boundary layer, and between the boundary layer and the free troposphere, can be exposed on the basis of mixing that follows a mixing-length hypothesis.However, when the boundary layer is heterogeneous or highly dynamic (such as during morning convective growth and evening transitions to stability, or when there are intermittent mixing processes) more detailed analysis is needed to evaluate the appropriateness of the mixing line approach.Figure 3 reveals periods when δD and CO 2 both feature sudden and large variations.During the night in particular, both tracers reveal the slow mixing is punctuated by bursts within the boundary layer or originating at the top of the observed profile.Such intermittent dynamical processes that affect the momentum and heat balance of the nocturnal boundary layer structure have been noticed in previous studies (Poulos et al., 2002;Salmond and McKendry, 2005).An example of the influence of these events on trace gases is seen at 10:00 p.m. on 15 February: δD and CO 2 suddenly decrease at 300 m and within the profile, suggesting a turbulent burst that leads to enhanced exchange with the free troposphere.Subsequently, δD decreases while CO 2 increases, suggesting a turbulent burst within the boundary layer that mixes depleted δD downward and transports the high CO 2 concentration upward, or from a lateral source.This type of transient dynamical feature is not described by a mixing line analysis.
The formation of a LLJ during the second night is a sustained dynamical structure which clearly illustrates limitations of a mixing line analysis.LLJs can result from many generation mechanisms.In this case, the very stable conditions, evident in both the tower data (Fig. 3a) and from the TLS (Fig. 7), tends to extinguish turbulence above the typically thin surface boundary layer, which leads to unbalanced horizontal pressure gradients that accelerates the flow until turbulence is reestablished (e.g., Businger, 1973).Figure 7 shows a series of very high resolution profiles during the time of the LLJ from the TLS.The turbulence profiles show the 10-15 m deep surface layer.Before the jet appears, the boundary layer profiles (the first seven profiles in Fig. 7) show a steep inversion in the lowest 15 m that indicates the depth of the surface layer, while the structure above is near neutral.The jet appears at 21:05 local time as a result of the enhanced stability.The gas tracers show that boundary layer air is replaced by well mixed air mass transported into the region above the surface layer by the jet, while air below the height of the jet retains CO 2 and δD values representative of the surface fluxes (Fig. 3d and e).From the wind speed data (Fig. 3b) the LLJ is identified mostly between 50 and 300 m.The relationship between the height of the LLJ and trace gas measurements from the Picarro analyser is partially masked by the inherent smoothing applied to the δD and H 2 O during the application of the memory correction (Appendix B).The three subsequent profiles (21:05, 21:20 and 21:30, Fig. 7) show the boundary layer becomes deeper allowing export of the surface air.The turbulent dissipation rate, ε, shows that a surface based turbulent boundary layer extends up to about 150 m.When the jet disappears, CO 2 and q are transported upward again at the slow rate governed by the nighttime stability.The advective influences of the LLJ on the moisture profile would invalidate the conditions required for the isotope GML method to be useful.The impact of a similar advective influence on the profile is illustrated well by the arrival of the frontal air mass in the evening of the 17 February (Fig. 8a).While reminiscent of the conceptual depiction in Fig. 1c, the vertical structure is not governed by a mixing process and so is not well modeled with a mixing line.Similarly, changes in the air mass over the depth of the profiles associated with the passage of the front on the 18th (Fig. 8b) are not captured by a mixing line analysis.Under these types of advective conditions the GML approach may be valid in the surface layer where the surface exchange dominates, but transient and dynamic phenomenon are not readily accounted for in a simple tracer mixing line analysis over different sections of the full 300 m profile.On the other hand, the isotopic information is useful in identifying these dynamic changes since they are easily identified where the data fail to conform to a simple mixing hypothesis.

Surface flux composition from mixing lines
Estimates of the δD of surface fluxes associated with evaporation and sublimation (i.e., δD F ) are obtained from the observed vertical profiles using a mixing line approach.The temporal mixing line (TML) analysis applied on temporal series in windows of 3.5 h of data taken from any single height yields noisy results and values of δD F that are not bounded by δD values of the snow and so appear unphysical.Testing different temporal windows from 1 to 6 h, applying smoothing filters and using different heights failed to produce results from the TLM method that could be deemed reliable  without a priori knowledge of the flux.This demonstrates that even over the time scale of a few hours, the dynamic and non-stationary nature of the boundary layer (illustrated in the previous section) prevents a simple two-box mixing model to capture the behavior of the profile.This result echoes similar concerns raised from continuous measurements elsewhere (e.g., Griffis et al., 2007).The shortcoming in the TML approach is because several fluxes affect the vapor δD at a given height during a particular time window (i.e., surface evaporation, mixing from above and regional horizontal advection frequently act at the same time).Figure 9a (green), shows flux values derived from the surface (10 m, and within the surface layer), which appear to be the most reliable because it is closest to the source and because the dynamics of the surface layer are simpler than the rest of the boundary layer.Uncertainty is shown in Fig. 9b as the standard deviation derived from the quadratic sum of precision of the measurement, uncertainty in all steps of the calibration and the estimate on the mixing line regression.The latter dominates.
Uncertainties in δD are associated with measurement pre- cision are ±5 ‰ at 1 g kg −1 and uncertainties arising from memory effects are typically ±2 ‰.This shows that while the TML method is not generally appropriate for boundary layer data, time series data obtained from short towers may be useful, especially over short time periods when the instantaneous gradient can be resolved (Good et al., 2012).The results suggests that studies using cryogenic methods to obtain samples that integrate over 10s of minutes to hours may be reliable when the conditions are chosen carefully.
The GML approach applied on profiles in the surface layer (taken from 0 to 60 m, Fig. 9 blue) yields values of δD F at noon and in the afternoon that remain between the predicted δD of snow sublimation and the predicted δD of evaporation of snow melt.It also captures the composition of the frost sample on the morning of 18 February.During other periods however, estimated δD F seems unphysical.This was most notable for the case of the frost event as the surface flux changes sign from upward during the day to downward at night in association with the frost.Similar limitations can be expected when the flux is small.The gradient mixing approach is sensitive to the height over which the mixing line is constructed.For instance, when calculated from 0 to 300 m, the approach yields significantly higher δD F values.The results show that while the GML method has advantages over the TML method it is not without limitations.The reason for some of the limitations and for the sensitivity to profile height warrants discussion.
The gradient mixing line approach requires a stationary profile.Most of the time, vertical profiles are not stationary, but rather reflect the transition between two (or more) states.surface fluxes lead to an enrichment of the vapor following a mixing line between night values and snow sublimation (Fig. 10a, red).The observed linear relationship suggests that the state is stationary, and explains the small sensitivity to profile height.This increases the confidence one may have in the δD F estimate.In contrast, later in the day, surface fluxes shift from sublimation to evaporation from standing surface liquid water, so that the surface end member varies with time.This transition leads to curvature of the mixing line in the evening (Fig. 10b, red).This explains why the δD F estimated from the full profiles are higher, and are probably larger than the true δD F value.Changes in the free tropospheric end member associated with shifts in synoptic-scale moisture origin similarly complicate the interpretation of the mixing line.
Figure 11 illustrates the behavior of mixing lines in the morning.During the night, the mixing line follows the same trajectory as that from the previous day and thus reflects the previous day's surface fluxes (Fig. 11a, black).When frost forms, however, the mixing line follows a new slope near the surface, which leads to curvature when the full profile is considered (Fig. 11a, red).This explains the particularly low estimates of δD F in the morning of 16, 17 and 18 February.These values do not represent the δD value of frost, but rather reflect the shift from one mixing line to another.When frost formation is stronger (relative to the mixing timescale) and dominates the profile, a new stationary mixing line is established and reflects frost formation very near the surface (Fig. 11b).Indeed the gradient mixing line approach applied to the bottom 60 m (i.e., the surface layer) can resolve the isotope ratio of frost early on 18 February.

Bounds on attribution of surface water flux
Surface snow δD is observed to increase by about 20 ‰ every day, and decrease during the night (Fig. 9).Since evaporative enrichment of pond water, enrichment during sublimation and re-equilibration of the snow pack with liquid give surface fluxes of different isotope ratio, the relevant processes can be constrained with the observed δD F value and the evolution of δD of the snow and melt ponds.The fluxes are decomposed by constraining the parameters of the mass balance model given in Sect.2.3.The model is initialized at 06:00 a.m. with the δD value of snow and water vapor set to the observed values of δD c = −175 ‰ and δD = −255 ‰ respectively.The observational constraints are typical δD values observed on 16 and 17 February at 06:00 p.m.: δD l = −140 ± 3 ‰ for ponds, δD c = −155±3 ‰ for snow and δD F = −220±15 ‰.Simulations were selected as plausible when the results were within one standard deviation of the observations.The mean of the ensemble of plausible solutions provides a measure of the most likely partitioning of fluxes, and the standard deviation and confidence intervals based on percentile bins provides an estimate of the uncertainty.Table 1 summarizes the results for each model parameter.
Simulation results show that both evaporative enrichment and re-equilibration with melt water are necessary to explain the observed approximate 20 ‰ enrichment of the snow each date.The role of snow re-equilibration is supported by observable changes in physical properties of snow over the course of the experiment from a "powder" texture to a larger grained and icier structure.This is consistent with significant melting and recrystallization.
Despite the wide range of possible values for the six model parameters, the set which produce simulations that match the observations depict a very consistent set of processes controlling the surface water budget (Fig. 12).Sublimation accounts for 70 to 73 % of the total surface water fluxes over the day, and for 59 to 60 % of the total snow loss.Snow samples contain 13 to 21 % of water that is refrozen liquid.The observed decrease of snow δD in the evening can be explained by drainage of evaporatively enriched liquid water that had accumulated within the snow pack during the day, but this is not resolved in the calculation.
A series of eight sensitivity experiments is used to test the importance of (1) the uncertainty in the observational constraints and (2) assumptions about the isotopic model on the resulting partitioning estimate.The sensitivity to observational uncertainty is tested by reducing the standard deviation of the constraints by a factor of two, and the importance of each constraint is assessed by withholding it from the calculation.Each test repeats the calculation with a 10 000 member ensemble to obtain optimal estimates.and values lower than one show that the uncertainty has decreased.
The sensitivity tests show that (1) the δD F constraint does not influence the uncertainty because it is itself uncertain, (2) the snow δD is the most significant constraint, and (3) when the pond water δD constraint is not used the uncertainty in the drainage fraction and pond evaporation is larger.The uncertainty of each of δD F , δD of snow and δD of pond water was reduced by a factor of two (tests 4, 5 and 6) and shows that there is only moderate change in the uncertainty on the final partitioning fractions which suggests that the measurement capabilities are sufficient.The importance of fractionation during sublimation is tested by allowing the fractionation during sublimation to vary between the model default value of 0 ‰ (e.g., Jouzel et al., 1987;Hoffmann et al., 1998) to 40 ‰ (the maximum value observed in laboratory experiments under windy conditions, Ekaykin et al., 2009).Adding this extra degree of freedom yields significantly larger uncertainty in the estimate of the partitioning.The degree of fractionation during sublimation is better constrained if the estimate of δD F is better, as illustrated by test 8 in which sublimation is allowed and the uncertainty of δD F is reduced by a factor of two.
The suite of sensitivity tests shows that the strongest constraint on the model of the snow pack is offered by measurements of the isotopic composition of the snow pack itself, ¾¹½¿% ¾¼¹ ¼% ¹ ¼% ¹¿¼% ½¹ % ½¼¼% Fig. 12.The partitioning of fluxes as described by an optimal fit of all model solutions to the observed isotopic constraints.Results are expressed as the percentage of the total surface water flux.The range is given as the maximum and minimum values from all simulations.The maximum likelihood estimates are given in Table 1.
and that uncertainty in the estimate of the partitioning between evaporation from liquid in the snow versus sublimation is much larger with only the estimate of the net flux.Because of this, the assumption on the existence of fractionation during sublimating (both equilibrium and kinetic) becomes an important aspect in obtaining robust estimates of the partition fractions.

Conclusions
Measurements of the lower boundary layer vertical profiles of specific humidity and the hydrogen isotope ratio of water vapor were made from a 300 m tall tower every 15 min during four winter days in Colorado.During the campaign we find that at the synoptic scale, the isotopic evolution reflects the origin and differing hydrologic histories of air masses.At the daily scale, the evolution of the isotope ratio reflects the sublimation of snow, evaporation of ponds and strong boundary layer mixing during the day.Since water vapor sources are largest and positive during the day and CO 2 emitted at the surface is trapped at night, they are complementary tracers of boundary layer mixing.Together, they show strong vertical mixing during the day and a shallow well stratified boundary layer during the night, in which mixing occurs mainly through intermittent bursts of turbulence.Several approaches were employed to deduce the isotope ratio of the flux using a mixing line analysis.Mixing line analyses applied to time series data yield very noisy results, due to the non-stationary and diverse set of processes (e.g., evaporation, turbulent mixing, slow shifts in tropospheric composition, etc.) taking place over the duration of the sections of data.Mixing lines constructed from vertical profiles yield more physically satisfying results during the day when thermal convection is well established, and the mixing line estimates from the surface layer (up to a few 10s of meters) provide the most reliable estimates, in agreement with other work (Williams et al., 2004;Griffis et al., 2007;Good et al., 2012).This confirms that measurements from short (10 m) towers are beneficial, but there remains a need for careful analysis of the gradient method if samples are taken sequentially rather than measuring the profile a given instant.The profile above the surface layer is influenced by more complex dynamics, and captures processes and source other than those associated with surface exchange, which negates the assumptions underlying the mixing model.The GML approach is also able to capture some frost formation events before sunrise.However, it fails during the night and every time the profiles deviate significantly from stationary state (e.g., when one of the end members shifts, or when the dominating process changes).Over the course of each day during the experiment, the progressive decrease in the δD values of the source vapor (i.e., δD F ) derived from the profiles suggests a shift from snow sublimation in the morning to pond evaporation through the day, which provides the evidence for physical changes in the characteristics of the source and a constraint for quantitative attribution.
Although limitations in mixing line analysis have been revealed at the processes level, the use of an optimal estimation approach allows synthesis of uncertain data from multiple sources with an adequate forward mass balance model to form an estimate of the flux partitioning.This is an important advance in the source attribution problem because it sets a path for more formal assimilation of surface isotopic data into detailed process models to constrain surface water and energy balance rather than relying on simple mixing methods.Additional information on water vapor source and exchange processes is provided by combining δD and δ 18 O, i.e. from deuterium excess measurements (e.g., Gat and Matsui, 1991).Developing appropriate calibration procedures for water vapor isotope ratio measurements remains a challenge, especially for deuterium excess when the humidity is low and highly variable as is typical in field settings.The analysis shows that while vapor measurements are important, measurements of the possible source reservoirs (especially soil water, snow pack, standing liquid and potentially precipitation, etc.) offer stronger constraints.These should be considered essential in developing observational programs that seek to use isotopic data for flux attribution.

Calibration of isotopic measurements
Raw δ (both D and 18 O) measurements are corrected to remove measurement dependence on mixing ratio and calibrated to the standard scale as shown schematically in Fig. A1.Calibration is based on data obtained from discrete injections of five known standards with a PAL autosampler and the Picarro vaporizer unit.Injections were made at mixing ratio values near 0.31, 0.62, 1.24, 3.10, 6.20 and 12.40 g kg −1 (using a 10 µm syringe for 6.20 and 12.4, and a 500 nl syringe for lower mixing ratios).Standard waters from Florida, Boulder, and the West Antarctic Ice Sheet, Greenland and Vostok were supplied by the University of Colorado Stable isotope Laboratory and tied to the SMOW-SLAP scale using IAEA primary standards.Humidity dependence is characterized by fitting a function, f , to measured values.Various functions can be used (geometric, splines, polynomials, etc.), and here we choose a cubic polynomial in the form δ = f (x) where x is the natural logarithm of q.The fit accounts for uncertainty in the δ and q using a Monte Carlo approach.The shape of the curve depends weakly on the δ value of the standard water, which is accounted for by linearly interpolating the regression coefficients between the known δ values of the standard water and the measured δvalue.The humidity correction is the difference between the δ-value at the measurement q-value and the δ-value at a reference value taken as q ref = 6.2 g kg −1 (Fig. A1a).Adjustment to the SMOW-SLAP scale is achieved with a quadratic fit between the values of standard waters known from mass spectrometer measurement and the measurements of those standard waters with the spectroscopic analyzer (Fig. A1b).A quadratic fit is needed to remove the non-linearity that would otherwise introduce errors of 2.8 ‰ for δD and 0.15 ‰ for δ 18 O.Uncertainty associated with instrument precision (which is a function of humidity), repeatability in measurement of standard waters, uncertainty in regression and curve fitting coefficients and accuracy of calibration waters is propagated though the approach using a Monte Carlo method to give an estimate on the measurement accuracy.The Monte Carlo ensemble is constructed by perturbing all quantities involved in the calibration by a normally distributed random amount proportional to one standard deviation of each quantity.The Monte Carlo ensemble has 10 000 members for each measured value, and the mean and standard deviation of the ensemble provides the estimate of the calibrated value and the accuracy on each measurement.The accuracy is from 6.6 to 3.3 ‰ for δD, and 1.0 to 0.4 ‰ for δ 18 O over the range of q from 1 to 4 g kg −1 .The one standard deviation precision of raw 0.16 Hz measurements is determined from laboratory measurements of water vapor with a fixed humidity and is from 5.5 and 2.6 ‰ for δD, and 1.0 to 0.30 ‰ for δ 18 O over the observed range of q.The accuracy of calibrated q measurements is 2 %.

Appendix B Algorithm for memory correction
Because of significant memory effects associated with low flow rates (30 cc min −1 ), the volume of the inlet and sample lines and instrument effects, the calibrated δ profiles for upward and downward motions of the elevator are different, and need to be reconciled by a posterior correction.Downward profiles have δ values typically 9.1 ± 9.5 ‰ lower than the upward profiles.This is because downward profiles retain the memory of the depleted vapor encountered at the top of the tower, whereas upward profiles retain the memory of enriched vapor encountered at the surface.
It is assumed that at each measurement time t, the composition of the vapor inside the isotopic analyzer, δ obs (t), is a combination of new sample vapor and the vapor from previous observations, and that the vapor is mixed with a time constant τ .The value of τ accounts for both the time required for mixing and the time scale of interactions between water molecules and all internal surfaces.Each vapor measurement reflects the environment vapor with a lag of t lag due to the transport time in the inlet.At each time t, δ obs is expressed as where δ env is the δ value of the environment and t is the time interval between measurements.Both τ and t lag vary with time, depend on wind, humidity and temperature conditions and on whether the elevator is ascending or descending.
The goal of the memory correction is to obtain an estimate of the time series of δ env (t).This estimation problem is regularized under the following assumptions: (1) τ (t) and t lag (t) are constant over each upward and downward motion of the elevator, (2) δ env varies slowly enough in time for δ env to be treated as constant for each of the 156 elevator cycles (a cycle includes one up and one down), and (3) the vertical structure of δ env (z) takes the form: This function is quite flexible, and allows a wide range of δ value profiles, with different average values, different vertical gradients, curvatures at the top or bottom, and up to two local maxima or minima.The time scale of the memory effect, of the order of several minutes, prevents statistically robust retrieval of any more detail in the profiles other than those included in Eq. (B2).This leaves 11 parameters to optimize for each elevator cycle: 7 parameters for the shape of δ env (z) and up and down values for τ and t lag .With 11 parameters over a 300 m profile it is reasonable to conceptualize an effective resolution of approximately 27 m, although the true resolution can be greater in regions of sharp gradients due to the structure of the assumed profile.Using all measurements from the δ env (t) time series for both the upward and downward motions of the elevator allows δ env (t) to be constrained efficiently.Minimizing the (root mean squared) difference between the upward and downward profiles effectively removes the memory effect.This mismatch is typically less than 2.5 ‰ for δD which is comparable in size to the measurement precision.

Fig. 1 .
Fig.1.Schematic depiction of three cases for applying mixing models to water vapor and δD profile measurements in the case of (a) a surface source with humidity anomalies that propagate upward, (b) a surface sink with humidity anomalies that originate at the bottom of the profile, or (c) exchange at the top of the profile.The tower is shown pictorially on the left with the instrument elevator moving either up or down.Curves on the left of each case indicate approximate water vapor or δD profiles, and boxes depict reservoirs that compose a mass balance.The quantities q and δD pertain to the measured profiles and q s and δD s are for the source or sink.Solid arrows show air mass exchange.Open arrows show water exchanges.The δD value of these fluxes is sought from a mixing line analysis.

Fig. 3 .
Fig. 3. Time height sections of (a) temperature ( • C), (b) wind speed (m s −1 ), (c) water vapor mixing ratio (i.e., q, g kg −1 ), (d) δD (‰), and (e) CO 2 volume mixing ratio (ppmv).Contours show potential temperature (K).Grey dots show where the Richardson number is the critical value, and the grey line estimates the depth of the convective layer.The time axis is in hours local Mountain Standard Time, and spans from 09:00 a.m. 15 February to 06:00 p.m. on 18 February.Arrows indicate examples of intermittent turbulent burst events.

Fig. 4 .
Fig. 4. Time series of (a) q(g kg −1 ), (b) δD (‰) and (c) CO 2 (ppmv) near the surface (taken at 10 m) and at the top of the tower over the four-day experiment.The time axis is in hours (local) Mountain Standard Time.Shading indicates the time between local sunset at 05:40 p.m. and sunrise at 06:50 a.m.

Fig. 5 .
Fig.5.Joint q − δD diagram showing tower data from 300 m (red and orange) and the base of the tower (blue and green).Data from the first three days of the experiment are shown as red and blue diamonds and the data from the final day of the experiment (after 07:00 p.m. MST on 17 February) are shown as orange and green squares.Theoretical curves are placed subjectively to provide context.The purple curve is a saturated Rayleigh curve assuming source vapor is initially at 70 % relative humidity over an ocean at 10 • C as might be typical of a North Pacific origin, and the cyan line shows a saturated Rayleigh curve assuming the source vapor is saturated over an ocean at 25 • C as is more typical of a Gulf of Mexico origin.The dashed cyan curve has the same origin, but the precipitation efficiency is 0.5.The orange line is for mixing between an air mass with q = 0.8 g kg −1 and δD = −350 ‰ and a source with δD = −190.

Fig. 6 .
Fig. 6.Profiles of δD for specific events that illustrate diurnal variations that can be captured by a mixing line analysis: (a) depletion associated with nocturnal mixing, (b) daytime enrichment from surface source, (c) frost formation.The red curve occurs several hours after the blue profile, and the evolution is indicated by the arrow.

Fig. 7 .
Fig. 7. Composite profiles from a high speed sensor on the tethered lift system during the night of 17 February 2010 during the formation of a low level jet.(a) The difference between potential temperature and potential temperature at the surface (K), (b) wind speed (m s −1 ) and (c) log of turbulence dissipation rate (m 2 s −2 ).Panels show a sequence approximately 10-30 min apart as indicated in (a) between 06:30 p.m. and 09:30 p.m. local time.

Fig. 8 .
Fig. 8. Profiles of δD for some events that demonstrate where a mixing line analysis fails: (a) a change in air mass origin associated with a change in the wind direction and, (b) decreases in the mean δD values of the entire profile associated with the passage of an air mass that has undergone precipitation.The red curve occurs several hours after the blue profile, as indicated by the arrow.

Fig. 9 .
Fig. 9. (a)Comparison between δD values of surface samples and estimates of δD F .The δD of snow samples (black circles) and values for vapor that would be in isotopic equilibrium with snow samples (red circles) provide upper and lower bounds for reasonable estimates of δD F .A blue asterisk shows the value of the measured frost sample and brown squares show mud samples.Estimates of δD F from gradient mixing line methods based on profiles from the surface up to 60 m (shown in blue) and from the surface to 300 m in magenta, and are usually within the expected range.Estimates of δD F from temporal mixing line approach from the data in 3.5 h windows from 10 m (green) are less robust.Only values where the regression is significant at the 95% confidence level are plotted.(b) Total error in the δD F estimate is shown as the standard deviation of the slope parameter in the linear regression, taking into account both total uncertainty in q and δD.Shading indicates time between sunset and sunrise.

Fig. 10 .
Fig. 10.Example of the evolution of profile mixing lines during the case of a surface water source.(a) Characteristic enrichment associated with sublimation of snow.(b) The case of a change in the isotopic composition of the source water from sublimation of snow during the morning to a source associated with a vapor flux from pond evaporation in the afternoon.Change in the background end member is associated with a slow enrichment of the higher altitude air mass.The black data points in (b) are the same as red data in (a).

Fig. 11 .
Fig. 11.As in Fig. 10 but for profile mixing line evolution during the formation of frost.(a) Light frost formation leads to a slight curvature in the data.(b) The case of strong frost shows the profile approaching the frost mixing line.The origin of high uncertainty in δD F is seen graphically as the sensitivity in estimating the intercept on the left axis from regression when the data is clustered on the right side of the diagram.
Fig. A1.Schematic depiction of calibration approach.(a) Raw measurements (δ raw ) are corrected for instrument dependence on humidity by applying an adjustment based on a (cubic) curve fit to standard water measured at different mixing ratio values.This standardizes δ measurements to a reference mixing ratio q ref = 6.2 g kg −1 .(b) Corrected values (δ cor ) are adjusted to the SMOW-SLAP scale to produce the calibrated values (δ cal ) using a calibration curve constructed as quadratic fit between values measured at q = q ref and known values of five working standards.Open symbols indicate measurements of standards and solid symbols show an observation.Solid arrows indicate the corrections.Values a and b are regression coefficients.All uncertainties are propagated using a Monte Carlo ensemble with 10 000 members for each measurement to give an estimate of accuracy.Error bars denote one standard deviation, and are greatly exaggerated to aid in the illustrative depiction.Dotted lines indicate the one standard deviation envelope on estimates of the calibration curves.Instrument precision is determined from laboratory tests.

Noone et al.: Water sources in the boundary layer
The importance of fractionation during sublimation is examined with sensitivity tests.Fourth, vertical heterogeneity in the snow pack could lead to different compositions uncovered with time.Frost formation on the existing snow pack would D.

Table 1 .
Table2reports results as the fractional change in the 25-75th percentiles range.Values shown in the table that are higher than one show where the possible range of values has increased, Optimal parameter estimates for partitioning in the surface water budget given in units of %.All fluxes are normalized by the total loss of snow mass.The sum f sub + f evapsnow + f evappond equals 100 % of the surface flux in each simulation in the Monte Carlo ensemble.Water budget fraction (%) f sub f evapsnow f evappond f melt f drain

Table 2 .
Fractional change in uncertainty of partition fractions from a series of eight sensitivity tests.Change in uncertainty is reported as the fractional change in the 25-75th percentile range.Values higher than one show that the range of values has increased, values lower than one means that the range of values is decreased.Tests 1, 2 and 3 remove each of δD F , δD of snow and δD of ponds.Test 4, 5 and 6 use uncertainty in those values reduced by a factor of two.Test 7 allows for fractionation during sublimation.Test 8 combines test 7 with reduced uncertainty in the estimate of δD F .