Effects of climatic seasonality on the isotopic composition of evaporating soil waters

Stable water isotopes are widely used in ecohydrology to trace the transport, storage, and mixing of water on its journey through landscapes and ecosystems. Evaporation leaves a characteristic signature on the isotopic composition of the water that is left behind, such that in dual-isotope space, evaporated waters plot below the local meteoric water line (LMWL) that characterizes precipitation. Soil and xylem water samples can often plot below the LMWL as well, suggesting that they have also been influenced by evaporation. These soil and xylem water samples frequently plot along linear trends in dual-isotope space. These trend lines are often termed “evaporation lines” and their intersection with the LMWL is often interpreted as the isotopic composition of the precipitation source water. Here we use numerical experiments based on established isotope fractionation theory to show that these trend lines are often by-products of the seasonality in evaporative fractionation and in the isotopic composition of precipitation. Thus, they are often not true evaporation lines, and, if interpreted as such, can yield highly biased estimates of the isotopic composition of the source water.

1 Introduction

Stable water isotopes (¹⁸O and ²H) are widely used in ecohydrology as tracers of the transport, storage, and mixing of water, from its origin as precipitation, through the soil, and ultimately to groundwater and streamflow (Kendall and McDonnell, 1998) or to plant uptake and transpiration (Dawson and Ehleringer, 1998). Water isotopes also reflect evaporation losses, through the progressive enrichment of ¹⁸O and ²H in the remaining liquid. Past applications of stable water isotopes in soil hydrology studies have focused on identifying evaporation fronts in the unsaturated zone (e.g., Allison and Barnes, 1983; Dawson and Ehleringer, 1998; Rothfuss et al., 2015; Sprenger et al., 2017a), quantifying groundwater recharge rates and mechanisms (e.g., Healy and Scanlon, 2010; McGuire and McDonnell, 2006; Sakakibara et al., 2017), quantifying root water uptake (e.g., Ehleringer and Dawson, 1992; Dawson et al., 2002; Volkmann et al., 2016; Rothfuss and Javaux, 2017), and identifying plant water sources (e.g., Brooks et al., 2010; Dawson and Ehleringer, 1991; Dawson and Simonin, 2011; Goldsmith et al., 2012; Evaristo et al., 2015; Hervé-Fernández et al., 2016; McCutcheon et al., 2017). A recent review by Sprenger et al. (2016) provides an extensive overview of isotope-based studies in the unsaturated zone. When expressed in the conventional δ notation and displayed together in a so-called dual-isotope plot (e.g., Fig. 1), variations in δ¹⁸O and δ²H in precipitation at any given location will typically exhibit a strong linear correlation (Craig, 1961; Dansgaard, 1964) termed the local meteoric water line (LMWL). As a water parcel evaporates, its isotopic composition will evolve along an evaporation line where the slope is determined by the relative evaporation rates of the different water isotopologues (Fig. 1). This evaporation line will generally have a shallower slope than the LMWL.

Soil and xylem water samples also often lie at an angle to the LMWL, and are frequently well described by linear fits (e.g., Brooks et al., 2010; Dawson and Simonin, 2011; Goldsmith et al., 2012; Evaristo et al., 2015). If these lines are evaporation lines, then extrapolating them to their intersection with the LMWL should yield the original composition of the pre-evaporation source water. Exactly this strategy has been used to infer source water compositions for soil water and xylem water (e.g., Evaristo et al., 2015; Hervé-Fernández et al., 2016; Javaux et al., 2016), as well as groundwater (e.g., Dogramaci et al., 2012) and stream water (e.g., Telmer and Veizer, 2000). Information on the source water composition is then typically used to draw conclusions about water cycling processes in the study systems. This inference should be valid if the evaporated samples all originate from a single source water.

But what if they don't? Is a linear trend, alone, sufficient evidence that the trend is actually an evaporation line? To date, no benchmark experiment has tested whether, and under what conditions, the trend line passing through fractionated soil water samples correctly identifies their source water.

Here we use simple numerical experiments, based on established isotope fractionation theory, to model the isotopic evolution of seasonally varying precipitation inputs, under the influence of seasonally varying evaporation processes. These simulations show that the resulting evaporated samples often fall along well-defined linear trends that are markedly different from evaporation lines, and therefore do not point to any meaningful source water composition.

2 Materials and methods

We simulated the isotopic composition of evaporating soil waters using equations based on the simple and widely used linear resistance model of Craig and Gordon (1965). We then introduced the effect of climatic seasonality by applying these equations to seasonally varying isotopic sources and atmospheric conditions.

2.1 Evaporative fractionation in open water and soils

The Craig and Gordon (1965) model estimates the joint effect of equilibrium and kinetic isotopic fractionation during the phase transition from liquid water to vapor. When resistance to transport in the liquid phase is neglected, the isotopic composition of the water vapor flux can be expressed as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{\mathit{\delta }}_{\mathrm{E}}={\displaystyle \frac{\left({\mathit{\delta }}_{\mathrm{L}}-{\mathit{\epsilon }}^{+}\right)/{\mathit{\alpha }}^{+}-h\cdot {\mathit{\delta }}_{\mathrm{A}}-{\mathit{\epsilon }}_{\mathrm{k}}}{\mathrm{1}-h+{\mathrm{10}}^{-\mathrm{3}}\cdot {\mathit{\epsilon }}_{\mathrm{k}}}},\end{array}$$

where δ_L and δ_A indicate the isotopic compositions of the evaporating surface and the atmosphere, h is the relative humidity of the atmosphere, α⁺ and ε⁺ are equilibrium fractionation factors, and ε_k is a kinetic fractionation factor. Here (and elsewhere in this paper), δ values and fractionation factors may refer to either hydrogen or oxygen isotopes unless otherwise noted. The δ notation expresses water isotope ratios as deviations, in parts per thousand, from Vienna Standard Mean Ocean Water (Kendall and Caldwell, 1998).

The equilibrium fractionation factor α⁺ (–) describes differences between the isotopic compositions of liquid and vapor phases at isotopic equilibrium, and is expressed here as the super-ratio of liquid to vapor isotope ratios. Its value is slightly larger than one, reflecting the fact that lighter molecules break their bonds more readily and thus are more abundant in the vapor phase. The values of α⁺ can be computed as a function of temperature T (K) using the well-established experimental results by Horita and Wesolowski (1994):

$$\begin{array}{ll}{\displaystyle }{\mathrm{10}}^{\mathrm{3}}\ln & {\displaystyle }\left[{\mathit{\alpha }}^{+}\left({}^{\mathrm{2}}\mathrm{H}\right)\right]=\mathrm{1158.8}\left(T^{\mathrm{3}}/{\mathrm{10}}^{\mathrm{9}}\right)-\mathrm{1620.1}\left(T^{\mathrm{2}}/{\mathrm{10}}^{\mathrm{6}}\right)\\ \text{(2)}& {\displaystyle }& {\displaystyle }+\mathrm{794.84}\left(T/{\mathrm{10}}^{\mathrm{3}}\right)-\mathrm{161.04}+\mathrm{2.9992}\left({\mathrm{10}}^{\mathrm{9}}/T^{\mathrm{3}}\right)\end{array}$$

$$\begin{array}{ll}{\displaystyle }{\mathrm{10}}^{\mathrm{3}}\ln & {\displaystyle }\left[{\mathit{\alpha }}^{+}\left({}^{\mathrm{18}}\mathrm{O}\right)\right]=-\mathrm{7.685}+\mathrm{6.7123}\left({\mathrm{10}}^{\mathrm{3}}/T\right)\\ \text{(3)}& {\displaystyle }& {\displaystyle }-\mathrm{1.6664}\left({\mathrm{10}}^{\mathrm{6}}/T^{\mathrm{2}}\right)+\mathrm{0.3504}\left({\mathrm{10}}^{\mathrm{9}}/T^{\mathrm{3}}\right).\end{array}$$

The equilibrium isotopic separation between liquid and vapor is then computed as ε⁺ = (α⁺ − 1) 10³ (‰).

The kinetic fractionation factor ε_k quantifies isotopic effects during net evaporation associated with the higher diffusivities of isotopically lighter molecules. Variations in ε_k are generally dominated by the relative humidity (h) of the air overlying the evaporating surface. Several expressions have been derived specifically for ε_k in soils (see Mathieu and Bariac, 1996; Soderberg et al., 2012). Here we use a simplified expression given as (Gat, 1996; Horita et al., 2008)

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }{\mathit{\epsilon }}_{\mathrm{k}}=\mathit{\theta }n(\mathrm{1}-h)\left(\mathrm{1}-D_{\mathrm{i}}/D\right){\mathrm{10}}^{\mathrm{3}}\left(\mathrm{‰}\right).\end{array}$$

The weighting term θ (–) accounts for the possible influence of the evaporation flux on the ambient moisture, and is usually assumed to equal one for small water bodies (Gat, 1996). The term D_i∕D is the ratio between the diffusivities of the heavy and light isotopes. Commonly accepted values are provided by Merlivat (1978): D_i∕D(²H) = 0.9755 and D_i∕D(¹⁸O) = 0.9723. The term n (–) accounts for the aerodynamic regime above the evaporating liquid–vapor interface. It ranges from n = 0.5 (fully turbulent transport that reduces kinetic fractionation, appropriate for lakes or saturated soil conditions) to n = 1 (fully diffusive transport, appropriate for very dry soil conditions). According to Eq. (2), in a dry atmosphere (h = 0), the kinetic fractionation factor is roughly 12.2–24.5 ‰ for ε_k(²H) and 13.8–27.7 ‰ for ε_k(¹⁸O).

We now consider the case of an isolated volume of water with initial isotopic composition δ₀ that evaporates into the atmosphere. As evaporation is the only flux, the remaining liquid volume decreases over time (a case sometimes referred to as a “desiccating” water body). We use x (–) to represent the fraction of the initial volume that has evaporated. The fraction remaining as liquid thus equals 1 − x. Assuming that the fractionation factors do not change during the evaporation process, the equation describing the isotopic composition of the residual liquid δ_L is (Gonfiantini, 1986)

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{\mathit{\delta }}_{\mathrm{L}}=\left({\mathit{\delta }}_{\mathrm{0}}-{\mathit{\delta }}^{*}\right)(\mathrm{1}-x{)}^m+{\mathit{\delta }}^{*},\end{array}$$

where δ^* [(‰) represents the limiting isotopic composition (i.e., the composition that the desiccating water volume would approach upon drying up) and the term m (–) is referred to as “temporal enrichment slope” (Gibson et al., 2016). These two terms can be computed as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }{\mathit{\delta }}^{*}=\left(h{\mathit{\delta }}_{\mathrm{A}}+{\mathit{\epsilon }}_{\mathrm{k}}+{\mathit{\epsilon }}^{+}/{\mathit{\alpha }}^{+}\right)/\left(h-{\mathrm{10}}^{-\mathrm{3}}\cdot \left({\mathit{\epsilon }}_{\mathrm{k}}+{\mathit{\epsilon }}^{+}/{\mathit{\alpha }}^{+}\right)\right)\end{array}$$

and

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }m=\left(h-{\mathrm{10}}^{-\mathrm{3}}\cdot \left({\mathit{\epsilon }}_{\mathrm{k}}+{\mathit{\epsilon }}^{+}/{\mathit{\alpha }}^{+}\right)\right)/\left(\mathrm{1}-h+{\mathrm{10}}^{-\mathrm{3}}\cdot {\mathit{\epsilon }}_{\mathrm{k}}\right).\end{array}$$

Equation (3) can represent an isolated volume of precipitation with initial isotopic composition δ_P that progressively evaporates into an atmosphere with an isotopic composition of δ_A. If the isotopic composition of the atmospheric vapor is unknown, it is common to assume that it is in equilibrium with precipitation (Gibson et al., 2008):

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{\mathit{\delta }}_{\mathrm{A}}=\left({\mathit{\delta }}_{\mathrm{P}}-{\mathit{\epsilon }}^{+}\right)/{\mathit{\alpha }}^{+}.\end{array}$$

As an introductory example, we modeled the isotopic evolution of an individual water volume by implementing Eqs. (1)–(6) with parameters T = 20 ^∘(C), h = 0.75 and n = 1. Figure 2 illustrates the resulting increase in the hydrogen and oxygen δ_L of the residual water during the evaporation process. The dual-isotope plot (Fig. 2c) shows the simultaneous behavior of the hydrogen and oxygen isotope ratios. As more of the water evaporates, the composition of the residual liquid gradually departs from the LMWL following a nearly linear trajectory. This trajectory is termed the evaporation line. Depending on the atmospheric parameters used in Eqs. (1)–(6), the slopes of evaporation lines will typically range from 2.5 to 5, markedly shallower than typical meteoric water lines, which usually have slopes of roughly 8 (Kendall and Caldwell, 1998).

Figure 2 Introductory example showing the evolution of the isotopic composition of residual water δ_L for the case of an isolated volume of precipitation that evaporates into the atmosphere. The initial composition (source water) δ₀ = δ_P is −6 ‰ (for δ¹⁸O) and −38 ‰ (for δ²H). (a) and (b) show the oxygen and hydrogen isotopic composition for increasing fractions of evaporation (decreasing fraction of residual liquid) as they approach the limiting composition, while (c) shows the same isotope effects in a dual isotope plot.

2.2 Accounting for the seasonality of atmospheric variables

The degree of evaporative fractionation will vary seasonally, reflecting seasonal changes in temperature and relative humidity. The isotopic composition of precipitation will also vary seasonally, reflecting seasonal shifts in moisture sources, air mass trajectories, and cloud processes (Craig and Gordon, 1965; Rozanski et al., 1993). With this in mind, we explore how these two seasonal patterns jointly shape the isotopic composition of the residual liquid remaining after rainfall partly evaporates from a soil.

We consider a 12 month period and for each month we use the mean isotopic composition of precipitation as source water for the model outlined above. Each month's precipitation then undergoes a seasonally varying amount of evaporation, and the isotopic composition of the residual water is determined separately for each month using Eqs. (2)–(6), along with that month's average temperature and relative humidity. In this approach, the isotopic composition of monthly residual water depends only on precipitation and evaporation during the same month. This simplified approach does not explicitly account for in-soil mixing processes, whose effects are discussed in Sect. 3.

We apply this approach to real-world weather and precipitation data from the Vienna Hohe Warte station, Austria. The full isotopic dataset is freely available, along with temperature and vapor pressure data, from the Global Network of Isotopes in Precipitation (GNIP Database), provided by IAEA/WMO and accessible at: https://nucleus.iaea.org/wiser. From GNIP we obtained long-term mean monthly values of precipitation δ¹⁸O at Hohe Warte, along with mean monthly air temperatures and vapor pressures, from which we calculated mean monthly relative humidities. Rather than using the precipitation δ²H values from GNIP, we instead computed them from δ¹⁸O using the equation for the LMWL at Hohe Warte: δ²H = 2.12 + 7.45δ¹⁸O. This ensured that each precipitation sample plotted exactly on the LMWL, thus aiding visualization.

The long-term mean monthly time series exhibit pronounced seasonality (Fig. 3). The seasonal temperature excursion is about 20 ^∘C, and monthly average δ¹⁸O ranges from −13 ‰ in winter to −6 ‰ in summer. The relative humidity ranges from roughly 0.85 in winter to 0.65 in spring and summer.

Figure 3 Long-term mean monthly air temperature (a), relative humidity (b) and oxygen isotopic composition in precipitation (c) for the station Vienna Hohe Warte, Austria.

To investigate the effect of seasonality in evaporation rates on residual liquid composition, we represented the evaporation-to-precipitation fraction (the variable x) by sinusoidal cycles with different amplitudes and timing. We did not consider transpiration fluxes, since the isotopic effects of transpiration are generally considered to be negligible. Moreover, to keep the example simple, we did not consider the seasonality of precipitation flux, although this could be easily included. The parameter n in Eq. (2) was fixed at 0.75 throughout the year. A numerical code to implement Eqs. (1)–(6) and apply them to the case of seasonal sources and climatic conditions is freely available at https://github.com/pbenettin/evaporation-lines and provided in the Supplement.

3 Results

3.1 Seasonal patterns in evaporated soil waters

The isotopic compositions of different source waters (mean monthly values of precipitation from the Vienna Hohe Warte station) and of the residual liquid water after evaporation (computed through Eqs. 2–6) are shown in dual-isotope space in Fig. 4. For this figure, we generated two hypothetical evaporation cycles, both peaking in July and having the same mean value $\stackrel{\mathrm{‾}}{x}$ of 0.10 (–) but with different degrees of seasonality. The weakly seasonal cycle had a peak-to-peak amplitude of just 0.02, and the more strongly seasonal cycle had an amplitude of 0.16. These x values are modest, representing conditions of limited evaporation that may be found in many temperate regions.

Figure 4 Effect of atmospheric seasonality on the isotopic composition of residual water from seasonally varying precipitation. The evaporatively fractionated residual water samples (green dots) cluster around a trend line (dashed line) which is much steeper than the individual evaporation lines (grey lines). The effects of strong and weak seasonality in evaporation rates (represented by x, the fraction of the initial volume that has evaporated) are shown in (a) and (b), respectively. The insets show the assumed annual cycles in evaporated fractions x.

The source waters (shown as yellow stars) vary along the LMWL reflecting the seasonal variability of atmospheric moisture sources and conditions, with isotopically lighter precipitation during colder months. The simulated residual water samples (shown as green dots) plot below the LMWL, with summer samples plotting farther from the LMWL than winter samples, reflecting their greater evaporative enrichment. The evaporation lines connecting individual source waters and residual waters are longer and shallower in summer than in winter, reflecting seasonal differences in temperature, relative humidity, and evaporated fraction x. As a result, the summer residual water samples plot farther away from the LMWL than the winter samples do, by an amount that reflects the seasonality in the evaporation process. The residual water samples follow a nearly linear trend (shown as a dashed line), which is markedly steeper than the evaporation lines for the individual source waters (shown as grey lines). The slopes of the evaporation lines range from 3.1 to 3.4; by contrast, the trend lines for the residual waters have slopes of 6.1 (Fig. 4a) and 7.1 (Fig. 4b), close to the assumed LMWL slope of 7.45. Note that whenever the residual water trend line has a slope that is close to that of the LMWL, the location of the intersection between these two lines will be highly uncertain.

Because the simulated residual water samples can be fitted easily with a simple trend line, it may seem logical to interpret this trend line as an evaporation line, and to infer an apparent source water end-member from its intersection with the LMWL. In the case of an isolated water parcel that is progressively evaporated (as in Fig. 2c), this approach could yield a reasonable estimate of the original source water. However, when residual water samples do not come from a single source, the trend line is not an evaporation line, and the intercept of this trend with the LMWL can lie far away from the average source water (Fig. 4a); the intercept can even lie far outside the range of all the source waters (Fig. 4b).

Figure 5 illustrates how different degrees of seasonality in evaporation patterns may yield different trend lines in residual water samples, with different intercepts with the LMWL. The individual source waters and evaporation lines are the same as in Fig. 4. The five trend lines in Fig. 5 are associated with different seasonal evaporation cycles, which feature similar low evaporation fractions in winter (roughly x = 0.04), but different evaporation fractions in summer (roughly x = 0.15 to x = 0.60). The evaporation cycles with higher summer peaks correspond to trend lines with shallower slopes and less negative intersections with the LMWL. All of the intersections lie far from the true mean source water; indeed none of them lie within the range of the individual monthly source waters.

Figure 5 Examples of trend lines and intercepts arising from various seasonal evaporation patterns (inset).

If the seasonal cycle of evaporative fractionation is not in phase with the seasonal cycle in source water composition (that is, if the most strongly fractionated sample is not also the one with the heaviest initial isotopic signature), the residual water samples will trace out a hysteresis loop. In Fig. 6, the source waters are the same as those in Fig. 4, but the seasonal evaporation cycle has been shifted by two months. The width of the resulting hysteresis loop depends on the amplitude of the seasonal cycle in evaporation, and how far out of phase it is with the seasonal cycle in precipitation isotopes. Even where such hysteresis loops exist in nature, they may be difficult to detect due to measurement uncertainties and environmental noise.

Figure 6 Hysteretic pattern arising in the computed residual liquid when the seasonality of evaporation rates is shifted relative to the seasonality of the source water composition (see inset).

3.2 Mixtures of evaporated soil waters

In Figs. 4–6, each residual water sample is derived from a discrete monthly precipitation source water sample. Real-world soil waters, by contrast, can be expected to contain mixtures of waters with different ages, and thus different source water signatures and evaporative fractionation trajectories. For simplicity, we simulated the soil as a well-mixed reservoir that integrates each month's residual waters. Mathematically this means that the composition of the soil pool is an exponentially weighted running average of the residual water samples shown in Fig. 4a. For purposes of illustration we used a time constant of six months, such that the same-month contribution to each sample is roughly 15 % and the contribution from the previous 12 months is roughly 86 % of the total. The results are shown in Fig. 7a.

Figure 7 Isotopic composition of modeled soil waters (triangles) obtained by mixing evaporated source waters through an exponential function. Soil water isotopic compositions were obtained starting from (a) long-term monthly sources (as in Figs. 4–6), and (b) individual monthly isotopic sources, recorded at Vienna Hohe Warte station between 1961 and 2015 (grey dots).

The same procedure was used to create Fig. 7b, except we considered each of the approximately 600 individual monthly δ¹⁸O and δ²H values available at Hohe Warte (the cloud of grey dots) as meteoric source waters. These source waters were not constrained to lie along the LMWL, in contrast to the analyses presented above. These source waters were individually evaporated and fractionated, by amounts that depended on the individual monthly temperature and relative humidity (and the same seasonal cycle in the evaporated fraction x that was assumed in Figs. 4a, 6, and 7a). We then applied the same running weighted time averaging used in Fig. 7a, with the resulting cloud of residual water samples shown in Fig. 7b. The more that the residual water samples are time-averaged, the more their scatter will be compressed and the smaller the portion of the dual-isotope plot they will occupy, but their trend line will remain almost the same. The exponentially weighted averaging used here also introduces a time lag of roughly 3–4 months between the seasonal cycle in the source water and the seasonal cycle in the time-averaged soil water. For this reason, the isotopically heaviest soil water samples are found in October even though the isotopically heaviest precipitation falls in the summer. (Different time constants in the weighted averaging would yield different lag intervals.) Similar lag periods of several months are often found in experimental studies (e.g., Sprenger et al., 2017b).

Due to the scatter among the source water samples in Fig. 7b, the evaporatively fractionated residual water samples are less collinear than in Fig. 7a. Nonetheless, in both cases the trend lines intersect the LMWL far from the true mean source water. Because the intersection point lies within the range of the individual winter precipitation samples, however, there is a risk that one could incorrectly infer that it represented a winter-precipitation source water for the evaporated soil samples (when in fact the winter precipitation in these simulations has hardly been evaporated at all).

4 Discussion and concluding remarks

The analyses presented above serve as a reminder that isotope samples need to be understood as combining the effects of source variation, mixing, and fractionation. Indeed, in our examples, all three of these effects jointly determine the isotopic patterns in the evaporated soil samples.

All else equal, the greater the isotopic variability in precipitation (and thus the larger the range of source waters), the closer the slope of the evaporated samples will lie to the LMWL (Fig. 8a). (Conversely, in the absence of any variability in precipitation, the evaporated samples would trace out evaporation lines instead.) All else equal, the greater the seasonality in evaporative fractionation, the more the slope of the evaporated samples will deviate from the LMWL (Fig. 8b). The intercept of the trend line with the LMWL is driven purely by these geometric considerations, and has no significance in its own right.

Figure 8 Effects of precipitation source variability (a) and evaporation variability (b) on the trend line that interpolates evaporated soil water samples (dots). The trend line is much steeper than the evaporation lines, unless there is no variability in the source water isotopic composition.

The seasonality of evaporative fractionation combines two factors: the variation in the slope of the evaporation line, and the variation in the amount of water lost to evaporation (as quantified by the evaporation fraction x), which determines how far out on the evaporation line the evaporated samples are found. In most real-world situations, the second of these factors is likely to have a greater influence on the trend line of the evaporated samples (and thus on its intersection with the LMWL).

The more mixing the evaporated samples undergo, the more their variability will be compressed. In addition, the more closely that the variations in precipitation isotopes and evaporation rates are synchronized, the more the evaporated samples will follow a trend line; conversely, if they are out of phase, they will form a hysteresis loop.

Our analysis of the effects of variability in source signatures and evaporative fractionation has been couched in terms of seasonal patterns, but similar considerations apply to variations at other timescales as well. For example, under more arid conditions the evaporation line will have a flatter slope and evaporative losses will be greater, both factors that will push evaporated samples farther from the LMWL. If those atmospheric conditions are also correlated with isotopically heavier source waters, the resulting residual water trend line will be similar to those we have simulated here (with a slope much steeper than a true evaporation line, and typically intersecting the LMWL far from the average source water).

The intersection between this trend line and the LMWL gives a heavily biased estimate of average source water, but what could give a better one? One can see from Fig. 7 that a reasonable estimate of the average source water could be obtained by translating the individual evaporated samples back to the LMWL along assumed evaporation lines, yielding estimates of their pre-evaporation compositions which are then averaged. Ideally the slope of each evaporation line would be determined from atmospheric conditions that are specific to each evaporated sample. But even where these are unknown, any reasonable estimate of the evaporation slope will yield much better results than the slope of the trend line through the evaporated samples.

We have chosen a relatively simple model to simulate the evaporative fractionation of the residual water samples. More sophisticated models of evaporative fractionation in soil water have been proposed (see, e.g., Mathieu and Bariac, 1996; Soderberg et al., 2012; Dubbert et al., 2013; Good et al., 2014). Results from these models may provide more accurate estimates of the kinetic fractionation factors and thus of the slope of the evaporation line (and its variability). However, use of these models is unlikely to yield qualitatively different results from those in Figs. 4–7, because other reasonable estimates of the fractionation factors (and thus of the slopes of the evaporation lines) will make little difference to the slope of the trend line running through the evaporated samples. Our analysis also invokes the simplifying assumption that (for example) July's rainfall only evaporates under July conditions. But if some of July's rainfall is stored until August, September, October, etc. then some of it should also evaporate under those conditions. Arguably our analysis could be superseded by a detailed process model that simulates the time-dependent storage and release of water in soils. However, such a model would complicate the analysis considerably and we have no reason to believe that it would yield substantially different results.

The data and equations presented in this paper are not novel, and many readers will not be surprised by our conclusion that trend lines through evaporated samples can differ widely from true evaporation lines. Nonetheless, our analysis shows how residual water trend lines can result from the interplay of seasonally varying isotopic inputs and evaporation rates, and shows that their intersection with the LMWL will generally be a highly unreliable guide to the average source water composition. Analyses that have used these trend lines to identify the compositions of source waters may be substantially in error, and therefore should be re-examined. Because plant uptake is generally not strongly fractionating, isotopic variations in soil water are likely to be transferred to plant xylem, and thus we expect that our conclusions will also apply to xylem water as well. Whether on seasonal or synoptic timescales, the regional and global energy dynamics that drive variations in source water composition and evaporative fractionation are likely to be widespread. Thus, although results for individual sites and time periods may differ in quantitative details from those presented here, we expect the qualitative patterns to be general.