Causes and consequences of pronounced variation in the isotope composition of plant xylem water

Stable isotopologues of water are widely used to derive relative root water uptake (RWU) profiles and average RWU depth in lignified plants. Uniform isotope composition of plant xylem water (δxyl) along the stem length of woody plants is a central assumption of the isotope tracing approach which has never been properly evaluated. Here we evaluate whether strong variation in δxyl within woody plants exists using empirical field observations from French Guiana, northwestern China, and Germany. In addition, supported by a mechanistic plant hydraulic model, we test hypotheses on how variation in δxyl can develop through the effects of diurnal variation in RWU, sap flux density, diffusion, and various other soil and plant parameters on the δxyl of woody plants. The hydrogen and oxygen isotope composition of plant xylem water shows strong temporal (i.e., sub-daily) and spatial (i.e., along the stem) variation ranging up to 25.2 ‰ and 6.8 ‰ for δ2H and δ18O, respectively, greatly exceeding the measurement error range in all evaluated datasets. Model explorations predict that significant δxyl variation could arise from diurnal RWU fluctuations and vertical soil water heterogeneity. Moreover, significant differences in δxyl emerge between individuals that differ only in sap flux densities or are monitored at different times or heights. This work shows a complex pattern of δxyl transport in the soil–root–xylem system which can be related to the dynamics of RWU by plants. These dynamics complicate the assessment of RWU when using stable water isotopologues but also open new opportunities to study drought responses to environmental drivers. We propose including the monitoring of sap flow and soil matric potential for more robust estimates of average RWU depth and expansion of attainable insights in plant drought strategies and responses. Published by Copernicus Publications on behalf of the European Geosciences Union. 4854 H. P. T. De Deurwaerder et al.: Variation in plant xylem water isotopic composition


Introduction
The use of the stable isotope composition of water has strengthened ecohydrology studies by providing insights into phenomena that are otherwise challenging to observe, such as relative root water uptake depth (RWU depth) , belowground water competition, and hydraulic lift (Hervé-Fernández et al., 2016;Meunier et al., 2017). Compared to root excavation, the technique is far less destructive and labor intensive. This makes it more flexible for studying multiple individuals across spatial and temporal scales (i.e., individual to ecosystem, daily to seasonal) (Dawson et al., 2002). Moreover, the study of stable isotope composition of xylem water measures the real effects of RWU at different depths, whereas excavation yields only root distribution and architecture. The advantages and wide applicability of this method make it a popular technique that pushes the boundaries of ecohydrology (Dawson et al., 2002;Lanning et al., 2020;Rothfuss and Javaux, 2017;Yang et al., 2010) A variety of methods are used to infer average RWU depth from the isotope composition of plant xylem water (δ xyl ), but all rely on a direct relationship between the isotopic compositions of plant xylem and soil water (Ehleringer and Dawson, 1992). All have two key assumptions. The first is that the isotope composition of plant xylem water remains unchanged during transport from root uptake to evaporative sites (e.g., leaves and non-lignified green branches). Hence, isotopic fractionation -i.e., processes that cause a shift in the relative abundances of the water isotopologues, driven by their differences in molecular mass -does not occur during the transport from the uptake to the evaporative site (Dawson et al., 2002;Dawson and Ehleringer, 1991;Walker and Richardson, 1991;Wershaw et al., 1966;White et al., 1985;Zhao et al., 2016;Zimmermann et al., 1967). Second, all methods assume that xylem water provides a well-mixed isotope composition of water from different soil layers; sampled xylem water instantaneously reflects the distribution and water uptake of the roots independent of the timing or height of sampling.
The first assumption is relatively well supported. Isotopic fractionation at root level does not raise concerns for most RWU assessments using water isotopologues  except for kinetic fractionation that might occur during water transport across the root membrane in extreme environments (Lin and Sternberg, 1993;Ellsworth and Williams, 2007;Zhao et al., 2016). Similarly, isotopic fractionation of water within an individual plant, although possible, is generally not considered a serious problem (Cernusak et al., 2005;Dawson and Ehleringer, 1993;Mamonov et al., 2007;Yakir, 1992;Zhao et al., 2016), but this perception was recently contested by Barbeta et al. (2020). However, the second assumption of time and space invariance of the isotope composition of xylem water has, to our knowledge, never been assessed.
Various plant physiological processes, ranging from very simple to more complex mechanisms, could have an influence within plant variation in δ xyl at short timescales, i.e., sub-daily to sub-hourly. For instance, plant transpiration during the day is regulated by stomata according to water supply and atmospheric demand and follows well-known diurnal patterns (Epila et al., 2017;Steppe and Lemeur, 2004). This results in a changing water potential gradient between soil and leaves throughout the day (Fig. 1a, b) which in turn affects the depth of the average RWU (Doussan et al., 2006;Goldstein et al., 1998;Huang et al., 2017). Hence, shifts in a plant's capacity to take up water at different soil layers during the day can generate diurnal variation in the mixture of isotope composition from water taken up from various depths (Fig. 1c). Subsequently, this water mixture moves up along the xylem with the velocity of the sap flux density. As these sap flux densities depend on species and individual-specific hydraulic traits and on their responses to atmospheric water demand and soil moisture availability, complex dynamics in isotopic composition will emerge and propagate through the plant. The above hypothesis, if true, would make the comparison of isotopic data among individuals, species, and studies difficult.
In this study, we provide a critical assessment of the assumption of δ xyl invariance along the length of woody plant stems and over short time periods. We first show that variation in δ xyl along the length of lignified plants exceeds the expected measurement error using three independent datasets including (i) canopy trees and lianas sampled at different heights in French Guiana, (ii) plant species in northwestern China (Zhao et al., 2014), and (iii) European beech and silver firs in southwest Germany (Magh et al., 2020). Second, we build a simple mechanistic model that incorporates basic plant hydraulic transport processes. The model predicts that diurnal changes in water potential gradient between soil and roots result in shifting sources of water absorption that differ in their isotope composition.  Figure 1. The use of stable water isotopes (δ 2 H and δ 18 O) to assess the relative depth of root water uptake (RWU) requires a depth gradient in isotopic composition of soil water (δ 2 H S ) to be present (a, sub-panel 1) as only then can the relative contribution of different soil layers to the isotopic composition in a plant's xylem water (δ 2 H X ) be derived. These δ 2 H S gradients occur naturally as the result of evaporative soil drying during drought conditions; however, these conditions also result in the formation of a gradient in soil matric potential ( S ), ensuring an increasing S with depth (a, sub-panel 2). RWU and sap flow in plants are passive processes in which water flows in the direction of decreasing water potentials. Specifically for RWU, this implies that water influx through the absorptive root area (A R ; a, sub-panel 3) of a plant's root is facilitated whenever the water potential in the root ( R ) is more negative than the surrounding S . As A R and S are generally not uniform with soil depth (z), the relative contribution of a specific soil layer to RWU will depend on (i) the difference between S and R in that soil layer and (ii) the relative amount of absorptive root area in that soil layer. Stable water isotope techniques assume that the δ 2 H X reflects the contribution of δ 2 H S from all soil layers. However, this does not account for the diurnal fluctuations in R which are invoked by the diurnal patterns in a plant's transpiratory water demands (b). Typically, more negative R values are observed when water demands are high, i.e., around midday. However, a decrease in R will result in higher RWU and alter the contribution of different soil layers to RWU. Specifically, dryer and shallower soil layers, with more negative S , could start contributing to RWU as R decreases (c). For example, in the early morning (situation t x ) when R is high, only deeper soil layers where S is greater than R contribute to the overall δ 2 H composition of the RWU flux. As L and R decrease towards midday (situation t y ), more water can be absorbed from shallower soil layers. As the A R in these shallow soil layers is high, they strongly affect the relative contribution of δ 2 H S entering the plant. Hence, diurnal fluctuations in R will result in fluctuating mixtures of δ 2 H S entering the plant. As these δ 2 H S mixtures are transported along the xylem pathway, they produce variances in δ 2 H X which could complicate RWU assessments via stable water isotope analysis. 0.15, 0.30, 0.45, 0.60, 0.90, 1.20, and 1.80 m) with a soil auger and in close vicinity to the sampled individuals. Samples were placed in glass collection vials, sealed with a cap, and frozen awaiting cryogenic vacuum distillation (CVD; 4 h at 105 • C). When the weight loss of a sample resulting from the extraction process was below 98 %, the sample was excluded (based on Araguás-Araguás et al., 1998) (see Fig. S1 in the Supplement).
The isotope composition of the water in the samples was measured with a wavelength-scanned cavity ring-down spectrometer (WS-CRDS; L2120-i, Picarro, California, USA) coupled with a vaporizing module (A0211 high-precision va-porizer) and a micro combustion module to avoid organic contamination (Martin-Gomez et al., 2015;Evaristo et al., 2016). Post-processing of raw δ readings into calibrated δ values (in ‰; V-SMOW) was performed using SICalib (version 2.16; Gröning, 2011). More details on the sampling site and sampling procedure can be found in Supplement methods A. Plant δ xyl was sampled at high temporal resolution in the Heihe River Basin (HRB), northwestern China. Four distinct study locations differing in altitude, climatological conditions, and ecosystem types were selected. At each location, the dominant tree, shrub, and/or herb species were considered for sampling. In August 2009, Populus euphratica was sampled in the Qidaoqiao riparian forest (42 • 01 N, 101 • 14 E) and Reaumuria soongarica in the Gobi Desert ecosystem (42 • 16 N, 101 • 17 E; 906-930 m a.s.l., above sea level). In June-September 2011, Picea crassifolia, Potentilla fruticosa, Polygonum viviparum, and Stipa capillata were measured in the Pailugou forest ecosystem (38 • 33 N, 100 • 18 E; 2700-2900 m a.s.l.). All species were sampled every 2 h over multiple days (3-4), except for P. crassifolia, which was measured hourly. Stem samples were collected for trees and shrubs, while root samples were obtained for the herb species. More details are available in Zhao et al. (2014).
Upon collection, all samples were placed in 8 mL collection bottles and frozen in the field stations before transportation to the laboratory for water extraction via CVD (Zhao et al., 2011). Both δ 18 O and δ 2 H were assessed with a Euro EA3000 element analyzer (Eurovector, Milan, Italy) coupled with an Isoprime isotope ratio mass spectrometer (Isoprime Ltd, UK) at the Heihe Key Laboratory of Ecohydrology and River Basin Science, Cold and Arid Regions Environmental and Engineering Research Institute. Internal laboratory references were used for calibration, resulting in a measurement precision of ±0.2 ‰ and ±1.0 ‰ for δ 18 O and δ 2 H, respectively.

Field data Germany: high temporal variation in δ xyl
A δ xyl monitoring campaign studying mature silver firs (Abies alba; n = 3) and European beeches (Fagus sylvatica; n = 3) was conducted during progressing drought conditions (6-11 July 2017) at the Freiamt field site in southwest Germany. Isotopic composition of xylem water was obtained from branch samples, which were collected every 2 h between 07:00 and 21:00 at the same height and canopy orientation in the sun crown. Branches were stripped of bark and phloem tissue. A Scholander pressure chamber (Scholander, 1966), which allowed concomitant registration of water potential of the sampled branches, was used to extract xylem water directly in the field (Rennenberg et al., 1996). Both δ 18 O and δ 2 H of branch samples were determined with a wavelength-scanned cavity ring-down spectrometer (L2130i, Picarro, California, USA), followed by data correction using ChemCorrect™ (Picarro, 2010). For more details, see Magh et al. (2020).

Field data normalization
To aid visual comparisons, we use normalized δ xyl values (β 2 H X and β 18 O X ) which describe the deviation of an individual sample from the average isotopic composition (a) along the height h of the stem or (b) over 1 day: where N is the number of sampled heights or time steps during 1 day.
2.2 Part B: model exploration

Model derivation
The expected δ xyl at different stem heights within a tree during the course of the day can be derived from plant and physical properties such as root length density, total fine root surface area, water potential gradients, and the isotope composition of soil water (Fig. 2). We call this the SWIFT model (i.e., Stable Water Isotopic Fluctuation within Trees). To derive the SWIFT model, we first describe the establishment of δ xyl entering the tree at the stem base via a multisource mixing model (Phillips and Gregg, 2003). We subsequently consider vertical water transport within the tree, which relates to the established sap flow pattern.
To ensure consistency and clarity in variable declarations, we maintain the following notation in the subscripts of variables: uppercase roman to distinguish the medium through which water travels (X for xylem, R for root, and S for soil) and lowercase for units of time and distance (h for stem height, t for time, and i for soil layer index). A comprehensive list of variables, definitions, and units is given in Table  1. A schematic representation of the model is provided in Fig. 2a. Note that the model presented here focuses on hydrogen isotopes (i.e., 2 H/ 1 H) but can easily be used to study oxygen isotopes (i.e., 18 O/ 16 O).

Isotope composition of plant xylem water at stem base
The δ 2 H composition of xylem water of an individual plant at stem base (δ 2 H X,0,t ) (i.e., height zero; h = 0 m; Fig. 2a) at time t can theoretically be derived by calculating a weighted average of water taken up from different soil depths (Phillips and Gregg, 2003). The root zone is divided into n discrete soil layers of equivalent thickness z. Here, we assume a constant δ 2 H composition of soil water (δ 2 H S,i ) over time in each soil layer, a reasonable assumption when isotopic measurements are conducted during rain-free periods, allowing the expression of δ 2 H X,0,t to be as follows: where f i,t is the fraction of water taken up at the ith soil layer ( Fig. 2a) defined as the following: where RWU i,t is the net amount of water entering and leaving the roots at time t in the ith soil layer (RWU i,t is defined positively when entering the root). The current representation of the model does not account for water loss via the root system nor for the mixing of the extracted water from different soil layers within the roots until the water enters the stem base. When tree capacitance is neglected, the sum of RWU i,t across the entire root zone is equal to the instantaneous sap flow at time t (SF t ): where k i is the plant-specific total soil-to-root conductance over soil layer i, X,0,t is the water potential (i.e., the hydraulic head) at the base of the plant stem, and S,i,t is the soil matric potential at the ith soil layer (Fig. 2a). Total plant water potential is generally defined as the sum of the solute, pressure, gravity, and matric potential. As long-distance water transport through the xylem is studied, the osmotic potential and the kinetic energy head can be assumed negligible (Früh and Kurth, 1999). The xylem pressure potential is represented by X,0,t . The term z i is the gravimetric water potential necessary to lift the water from depth z i to the base of the stem, assuming a hydrostatic gradient in the transporting roots. The model considers z i to be a positive value (zero at the surface); thus z i is subtracted from S,i,t . A R,i is the absorptive root area distribution over soil layer i (Fig. 2a). This parameter, A R,i , can be derived from plant allometric relations with stem diameter (Čermák et al., 2006) and subsequently distributed over the different soil layers, considering the power-law distribution of Jackson et al. (1995). The total soil-to-root conductance is calculated assuming that the root and soil resistances are connected in series ( Fig. 2a): where k R is the effective root radial conductivity (assumed constant and uniform) and k S = K S,i / l is the conductance associated with the radial water flow between soil and root surface. The equation l = 0.53/ √ π × B i represents the effective radial pathway length of water flow between bulk soil and root surface (De Jong van Lier et al., 2008;Vogel et al., 2013) with B i giving the overall root length density distribution per unit of soil. K S,i is the soil hydraulic conductivity for each soil depth. K S,i depends on soil water moisture and thus relates to the soil matric potential S,i,t of the soil layer where the water is extracted. K S,i is computed using the Clapp and Hornberger (1978) formulation: where K s,max is the soil conductivity at saturation and b and sat are empirical constants that depend on soil type (here considered constant over all soil layers).
Subsequently, f i,t can be restructured as follows: where the root water to soil matric potential gradient is represented as (2) and (7) then allows the derivation of δ 2 H X,0,t as follows: This equation requires estimates of i,t , which is preferably measured instantaneously in the field (i.e., via stem and soil psychrometers for X,0,t and S,i,t respectively). However, as measurements of X,0,t are not always available, es-timatedˆ X,0,t can be derived from sap flow by reorganizing Eq. (4) intô which then allows the replacement of X,0,t withˆ X,0,t in Eq. (8).

Height-dependent isotope composition of plant xylem water
In our model, the water isotopologues simply move upwards from the stem base with the sap flow velocity. Assuming negligible diffusion, the δ 2 H isotope composition in xylem water at height h and time t (δ 2 H X,h,t ) is then the isotope composition of xylem water at stem base at time t − τ : where τ is the lag before δ 2 H X,0,t reaches stem height h (Fig. 2a), which depends on the true sap flux density in the xylem (SF V ). True sap flux density indicates the real speed of vertical water displacement within a plant and is derived by dividing SF t over the lumen area of the plant (A x ; Fig. 2a) i.e., the total cross-sectional area of the vessels. The τ value can be obtained from the mass conservation equality: Note that since most scientific studies express sap flux density as the sap flow over the total sapwood area (SF S ) rather than over the total vessel lumen area (SF V ), for consistency, we will present the model outputs as functions of SF S . Note that SF V presents the sap flux density normalized over the total vessel lumen area, and, as vessel lumen area correlates with plant diameter at breast height (DBH), there is no need for explicit consideration of DBH in the model for comparison among field measurements.
Model analyses show that the impact of the mutual diffusion coefficient of heavy water in normal water on the transport flux is negligible for plants with high sap flux densities, which is the case for the theoretical examples below. However, in plants with low sap flow densities, consideration of diffusion might be required. Diffusion might also be generated by water passing through a complex network of vessels, analogous to diffusion in a porous medium (see Supplement methods B for some analytical results and simulated cases of and a detailed discussion on the role of diffusion). SWIFT was implemented in R version 3.4.0 (R Core Team, 2017) and is publicly available (see GitHub repository HannesD-eDeurwaerder/SWIFT).

Model parameterization and analyses
The model's primary purpose is to gain insight into (1) which processes are capable of generating δ xyl variance and (2) how sensitive the variance in δ xyl along the stem is in response to the modeled plant hydraulic processes. To this end, we adopted the basic plant parameters from Huang et al. (2017) who studied soil-plant hydrodynamics of loblolly pine (Pinus taeda L.) during a 30 day extended dry down period (Table S1 in the Supplement). We started with synthetic basal sap flow patterns and volumes extracted from the model runs of Huang et al. (2017) for a typical drought day (day 11). Both basal sap flow patterns and volumes are repeated over the studied period as no variation between days is assumed. Sap flow follows the plant's water demand which is the result of daily cycles of transpiration driven by photosynthetic active solar radiation (PAR), vapor pressure deficit (VPD), and optimal stomatal response (Epila et al., 2017). Secondly, both the soil matric potential ( S,i,t ) and δ 2 H composition of soil water (δ 2 H S,i ) profiles with soil depth were adopted from Meißner et al. (2012) (Fig. S8; see Table S1 for equations) as driver data of the model and were assumed to stay constant over time. Since the measurements of Meißner et al. (2012) were conducted at a silt loam plot in the temperate climate of central Germany, corresponding soil parameters were selected from Clapp and Hornberger (1978). Subsequently, the following model simulations were executed (see Fig. 2a).
2. Analysis A2 studies the impact of temporal SF t variation at different tree heights and temporal patterns in δ 2 H X within a tree at various sampling heights (5, 10, and 15 m).
3. Analysis A3 studies the impact of temporal SF t variation on the isotope composition of xylem water and the timing of sampling. The representation of the profile of δ 2 H X along the full height of a tree is measured at different sampling times (09:00 and 11:00) with the standard parameterization given in Table S1.
4. Analysis B studies the variation in δ 2 H X due to differences in absolute daily average sap flow speed and diurnal patterns in the δ 2 H X in trees that differ solely in daily averaged SF V , which is set to 0.64, 0.42, and 0.19 m h −1 corresponding to SF S values of 0.09, 0.06, and 0.03 m h −1 , respectively.
All parameters of the four analyses are given in Table S1. The model simulations for each analysis were compared with a null model.

The null model
The null model adopts the standard assumption of zero variation in δ xyl along the length of the plant body but allows for potential measurement errors related to the extraction protocol. In reality, empirically obtained data will have some variation as observed values (Obs.δ xyl ) are the sum of the true δ xyl values and their extraction error (error_extraction): Obs. δ xyl = True δ xyl + error extraction .
Hence, the null model attributes any variance in isotopic composition to extraction errors with maximum extraction error ranges of 3 ‰ for δ 2 H samples (0.3 ‰ for δ 18 O) expected for water extraction recovery rates higher than 98 % (e.g., Orlowski et al., 2013). These extraction errors are negatively skewed following the Rayleigh distillation model which predicts that extraction error for incomplete water recovery will be negative and therefore Obs. δ xyl ≤ True δ xyl . The null model represents this error_extraction by a negative skew normal distribution (with location parameter ξ = 0 ‰, the scale ω = 3 ‰ for δ 2 H or 0.3 ‰ for δ 18 O, and shape α = −∞) (Azzalini, 2013).

Estimation of average RWU depth
Average RWU depths (i.e., the weighted mean of the depths of RWU with the uptake fractions at the different depths as weights) were derived from the simulated δ 2 H X values through the use of both the direct inference method and the end-member mixing analysis method. Together, these techniques represent 96 % of the applied methods in the literature , and the reader is referred to Rothfuss and Javaux (2017) for a complete discussion of both techniques. In line with the general approach assessing RWU with stable water isotopes, the average RWU depth is obtained by relating the δ 2 H X with the δ 2 H S,i depth profile. We compared average RWU depth estimates obtained from simulated δ 2 H X , as described in the analyses above, with the true average RWU depth. Here, the true average RWU depth was defined as the depth corresponding to the daily weighted average δ 2 H X and calculated as the weighted sum of δ 2 H X,i,t and the relative fraction of water taken up at each depth.

Transport dynamics and sensitivity analysis
We perform a basic model validation of our model assumption that the propagation of an isotopic signature is driven by diurnal sap flow dynamics and diffusion alone. In essence, the model assumes that once water with a given isotopic signature enters the stem, it moves upwards with the speed of sap flow and changes only due to the effect of diffusion.
The effects of capacitance on δ 2 H X dynamics by the release of storage water in the xylem flow can be ignored. To validate this assumption we compare model predictions against observed δ 2 H X dynamics monitored within a pine tree (Pinus pinea L.) following 2 H-enrichment in a controlled greenhouse experiment, as detailed in Marshall et al. (2020). δ 2 H X was measured at two heights (0.15 and 0.65 m) using a novel in situ technique, the borehole equilibration method. The model simulations performed consider the absolute ranges of sap flux densities during the entire monitoring campaign with the account of tree tapering effects on sap flux densities over the studied stem length (Supplement method C). Validation of diurnal variation in δ 2 H X requires high temporal resolution monitoring of δ 2 H X dynamics in plant stems with simultaneous high temporal resolution monitoring and characterization of sap flow, soil water potential, and isotopic composition. Such data do not yet exist to the best of our knowledge.
In addition, we performed two sensitivity analyses to assess the relative importance of each parameter in generating variance in δ 2 H X along the length of a plant. In both sensitivity analyses, we varied model parameters one at a time to assess the local sensitivity of the model outputs for soil type, sap flux density, root properties, and sampling strategies. The sensitivity analysis provides insight into possibilities for improving the design of field protocols by revealing potential key measurements and caveats in field setups. More details on the performed sensitivity analysis and validation of transport dynamics are available in Supplement method C.

Part A: empirical exploration
The null model assumes the constant isotopic composition of root water uptake with only limited variance in isotopic composition introduced by extraction errors (β 2 H X < 3 ‰; δ 18 O X < 0.3 ‰). However, pronounced δ 2 H X variance within individual plants exceeding the null model ranges, are observed in all three independent datasets. The normalized δ 2 H composition in xylem water (β 2 H X ) along the stem length of lianas and trees in French Guiana exceeded the null model by a factor of 3.2 and 4.3, respectively (Figs. 3c, S2). Differences of up to 13.1 ‰ and 18.3 ‰ in δ 2 H and 1.3 ‰ and 2.2 ‰ in δ 18 O were observed in individuals of trees and lianas, respectively (Supplement method A, Table A).

Part B: model exploration
Isotope composition of xylem water at stem base and basic model behavior At the stem base, simulated δ 2 H X,0,t displays a diurnal fluctuation (Figs. 2b, S7) that corresponds to the daily sap flow pattern (Fig. S7). This pattern is caused by a shifting diurnal average RWU depth. Early in the morning, when transpiration is low, most of the RWU occurs in deeper layers where soil matric potential is less negative and where soil water is more depleted in 2 H compared to the soil layers above ( Fig. S8a-b). As transpiration increases during the day, a significant proportion of RWU can now be extracted from the drier, shallower layers where the soil water is enriched in 2 H, i.e." having a higher δ 2 H. In the afternoon, as transpiration declines, the isotopic composition reflects again the composition of the more depleted soil water in the deeper soil layers, and it remains constant throughout the night because apart from diffusion, SWIFT does not consider mixing of the internal stem water. The mixing effects of diffusion are only noticeable at low sap flow speeds (Fig. 3b).
The highest δ 2 H X values (approx. −59 ‰) are found in alignment with the diurnal minimum of X,0,t (approx. −0.85 MPa; Fig. S7). At this moment, the difference between X,0,t and S,i,t is at a maximum, enabling water extraction from the upper and driest soil layers. Most root biomass is located near the surface (cf. Jackson et al., 1995;Fig. S8c), and uptake in these layers will result in relatively high contributions to the total RWU.
In contrast, differences between X,0,t and S,i,t are smaller in the early morning and late afternoon causing root water uptake in the upper soil layers to halt. The decreasing absolute range of i,t translates into higher proportions of RWU originating from deeper, more depleted soil layers. This causes δ 2 H X to drop to a baseline of approx. −67 ‰. This afternoon depletion of δ 2 H X will henceforth be referred to as the δ 2 H X -baseline drop.
Isotope composition of xylem water at different times, heights, and SF V Temporal fluctuation in δ 2 H X within a tree at 1.3 m (i.e., the standard sampling height; analysis A1; Fig. 2a) and other potential sampling heights (e.g., branch collection; analysis A2;   Table S1. (c) Field measurements of β 2 H X for six lianas ( ) and six trees ( ). Error whiskers are the combination of potential extraction and measurement errors of the isotope analyzer. A species-specific breakdown of the field data is provided in Fig. S2. The horizontal gray-colored envelope in all panels delineates the acceptable variance from the stem mean according to the null model (H 0 ), i.e., assuming no variance along the length of a lignified plant aside from the potential extraction error (i.e., 3 ‰). Herein, the dark gray envelope indicates the confidence interval comprising 95 % of potential extraction error (CI95). ment and the corresponding time needed to move the water along the xylem conduits. Note that the selected temporal resolution depends on whether the δ 2 H X -baseline drop at a given height equals the (stem base) minimum (here 1 min; see Fig. S12). In addition to sampling height, analysis A3 depicts the importance of sampling time (Fig. 3a). Outputs of analysis B predict that the occurrence and width of the δ 2 H X -baseline drop are a function of the sap flow velocity SF V (Fig. 3b). To aid model interpretation and comparability with field data, we (i) provide an illustrative example of normalized δ 2 H isotope composition of model-simulated xylem water (β 2 H X ) with consideration of extraction error (Fig. 4a) and (ii) display the relation between δ 2 H X variance and cumulative sap flow volumes, for which the piston flow dynamics in SWIFT originate from the lateral translation of the δ 2 H X fluctuation at δ 2 H X,0,t (Fig. 2b).

Potential biases in average RWU depth estimation
Both the timing of measurement (Fig. 5a) and SF V (Fig. 5b) influence average RWU depth estimates derived via the direct inference and end-member mixing analysis method (Fig. S9). Collection of tree samples at 1.30 m can result in erroneous estimation, deviating by up to 104 % from the average daily RWU depth (Fig. 5). Plotting the relative error in average RWU depth as a function of time and SF V (Fig. 5) shows that it is possible to time δ 2 H X measurements in a fashion that captures unbiased estimates of the average RWU depth. Xylem water sampling should be timed to capture the δ 2 H X that corresponds to water extracted at peak RWU, and the expected sampling time can be derived by considering the time needed for the water to reach the point of measurement (i.e., at 1.30 m in Fig. 5).

Transport dynamics and sensitivity analysis
Our sensitivity analyses show that the expected absolute error in average RWU depth assessment is directly related to both (1) maximum variance in and (2) the probability of sampling nonrepresentative δ 2 H X values. The maximum variance depends on the height, while the probability of sampling nonrepresentative areas depends on the width of the δ 2 H X -baseline drop (defined above). Hence, variation in δ 2 H X is determined by several factors, including the sampling strategy (timing and height of sampling), sap flow velocity (Fig. S10), and the belowground biophysical parameter (Fig. S11). We summarized the most important variables and 5 m (green) sampling heights with the formula provided. Thicker lines indicate model simulations without error, and lines connected by dots indicate a scenario of hourly sampling with consideration of the extraction error (i.e., a negative skew normal distribution; ξ = 0 ‰, the scale ω = 3 ‰, and shape α = −∞). (b-e) High temporal field measurements of β 2 H X of (b) two shrubs, (c) two trees, and (d) two herb species sampled in the Heihe River Basin (northwestern China) and (e) two tree species sampled in the Freiamt field site in southwest Germany. The horizontal gray-colored envelope in all panels delineates the acceptable variance from the stem mean according to the null model (H 0 ), i.e., assuming no variance along the length of a lignified plant aside from the potential extraction error (i.e., 3 ‰). Herein, the dark gray envelope indicates the confidence interval comprising 95 % of potential extraction error (CI95). A breakdown of the field data on species and individual level is provided in the Supplement figures (Figs. S3-S6).
as predicted by SWIFT which should be considered in subsequent RWU studies.
Plants on loamy soils show larger diurnal δ 2 H X variances in comparison to those on clay soils for a similar prevailing isotope gradient across the soil profile. Larger variances correspond to potentially larger errors, but the steeper slope of the δ 2 H X curve results in a thinner δ 2 H X -baseline drop. Hence, loamy soil can result in the potentially large error, but this is mediated by a lower probability of sampling nonrepresentative δ 2 H X values during the day.
The volume of water taken up by the plant (SF t ; Fig. S11b) affects xylem water potential of the plant at stem base (ˆ X,0,t ). Higher SF t requires more negativeˆ X,0,t , enabling the plant to access the enriched soil water of more shallow soil layers. Therefore, an increase in SF t results in the increase in maximum δ 2 H X values (increased maximum error) but also results in a smaller width of the baseline drop ( Fig. 2-3). Lower SF t results in smaller errors but a larger probability of sampling a nonrepresentative area (Fig. 3b).
Root properties, i.e., root membrane permeability (Fig. S11c), strongly influence both the total range of δ 2 H X variance and the width of the δ 2 H X -baseline drops. Decreasing root membrane permeability but with no alterations to the sap flow volumes results in thinner δ 2 H X -baseline drops but higher maximum δ 2 H X variance.
In addition, the true sap flow velocity (SF t per unit of lumen area) will determine the relative importance of diffusion on the δ 2 H X dynamics. Diffusion can cause a smoothing of the peak and a consequent increase in the width of the δ 2 H Xbaseline drop. However, as diffusion is proportional to the time that water isotopologues remain in the xylem, its absolute impact on δ 2 H X is negligible in plants with a high true sap flow velocity. In contrast, the impact of diffusion on δ 2 H X dynamics is substantial for plants with very low velocities, where water takes many days to pass from roots to leaves (see Supplement method B).
The role of diffusion was investigated using a stepwise 2 H enrichment experiment in Marshall et al. (2020) (Fig. 6). Analytical solutions of an advection-diffusion equation show that at 0.15 cm, a relatively small diffusivity was required to reproduce the initial increase in xylem isotope signature with values comparable to these reported for diffusivity of heavy water (Meng et al., 2018). However, at 65 cm, the value of diffusivity required to match the observed initial increase was much higher, suggesting other processes besides molecular diffusivity might contribute to the isotope transport (e.g., variable flow velocities within vessels and among vessels of the xylem network). Note also that the analytical solutions were not able to recover the second part of the curve where the isotope reaches the asymptotic enriched value, which is more gradual in the observations (Fig. 6). This also suggests a complex transport of δ 2 H X in the xylem.

Dynamic diurnal isotope compositions of xylem water along plant stems
Empirical field data show pronounced δ xyl variance along the stem length (Fig. 3) and over a sub-daily time period (Fig. 4).
Our model explorations suggest that basic plant hydraulic functioning can result in shifting mixtures of δ 2 H X entering the plant (Figs. 2-3). Daily X,0,t fluctuations interact with the S,i,t profile causing different parts of the root distribution to be active during the day. The fluctuations in δ 2 H X at the stem base propagate along the xylem with a velocity proportional to the sap flow, and this produces variability in sampled δ 2 H X that is much larger than the expected measuring error. Consequently, rather than being static, δ 2 H X values along the height of a plant should be envisioned as a dynamic diurnal process. Importantly, we show that high variance in δ 2 H X can result in an incorrect assessment of differences in average RWU depths between plants. Differences do not necessarily result from variability in average RWU depth but may result from monitoring plants at different heights (Figs. 2-3) and at dif-ferent times (Fig. 3a) or by comparing individuals which have different SF V (Fig. 3b) and xylem anatomical properties. For example, depending on SF V and the lumen area, the isotopic signal can take hours or days to travel from roots to leaves -as was also observed experimentally (Magh et al., 2020;Marshall et al., 2020;Steppe et al., 2010).
Low SF V allows multiple δ 2 H X -baseline drops over the length of a single tree. Sampled δ 2 H X can reflect the soil isotopic composition of the past several days. Our sensitivity analysis reveals that various soil and plant characteristics have an important role in determining both the daily maximum δ 2 H X variance, as well as the width of the δ 2 H X -baseline drop. These two characteristics directly impact (i) the expected maximum bias in estimates of average RWU depth and (ii) the chance of measuring δ 2 H X values that do not represent a mixture of all rooting layers during peak RWU (i.e., measurements in the baseline drop). Ultimately, these factors will challenge the use of stable water isotopes to study the terrestrial water fluxes, as recently reviewed by Penna et al. (2018). We additionally advocate that future research should explore the minimum set of (bio)physiological drivers and processes that require quantification to correctly interpret δ 2 H X along the hydraulic pathway length of a plant. Figure 6. Basic model validation comparing continuous in situ δ 2 H X measurements of a stepwise 2 H enrichment experiment (Marshall et al., 2020) with analytical solutions of an advection-diffusion equation at heights 0.15 m ( ) and 0.65 m ( ) on a pine tree (Pinus pinea L.). The source water of the intact-root, isotopic-enrichment greenhouse experiment is presented in gray. Model parameters, velocity, and diffusion were fitted by visual inspection independently for the two heights to match the initial increase in isotope signature (values reported in the bottom right).

General applicability of model and results
A necessary condition for diurnal shifts in RWU is the existence of water potential differences, e.g., more negative water potentials in the upper layers where trees usually have higher root density, which can cause a disproportional partitioning of diurnal RWU between deep and shallow roots over a diurnal course. The pronounced variance in δ xyl identified in this study is intrinsic to the isotopic tracing technique for RWU assessment as this method relies on the existence of a soil water isotopic profile. Such profiles are the result of soil evaporation, a process inextricably coupled to water potential heterogeneity, and hence to variance in δ 2 H X .
Plant transpiration results from a complex interaction between atmospheric demands (i.e., driven by VPD and radiation) and stomatal conductance that depends on tolerance for drought stress and soil moisture content. We may expect that diurnal fluctuation in radiation and VPD, and hence in water transport and depth of water absorption as modeled here, is a general phenomenon in nature. Moreover, much greater fluctuations in VPD and radiation should be expected under more natural conditions than the diurnal cycle described here, and these will increase the variability of transpiration fluxes, leading to even more complex dynamics of X,0,t . Specifically, the model simulations suggest that intraindividual variability of δ 2 H X will reflect the past changes in RWU dynamics, including RWU dynamics driven by changes in environmental demands. For instance, a changing degree in cloud cover that impacts sap flow dynamics can influence X,0,t rather abruptly (e.g., in lianas; Chen et al., 2015) and lead to instantaneous changes in the δ 2 H composition of the water mixture taken up at the root level. This can compli-cate the comparison of different plants sampled at different heights and times.
Note that, based on our model, we expect that soil isotopic enrichment experiments will generate extensive δ 2 H X variation along the length of trees whenever diurnal RWU fluctuations cause water extraction to shift between labeled and unlabeled soil layers. Furthermore, when enrichment experiments target trees with different hydraulic properties (such as SF V ) care should be taken to determine when and where to sample these trees to assess a potential 2 H enrichment in xylem water (Fig. 6, but see Magh et al., 2020;).

Alternative causes of δ xyl fluctuation
The SWIFT model provides a simple traceable and mechanistic explanation, using diurnal variations in SF t and RWU, for the pronounced variance and dynamic nature of the δ xyl fluctuations with plant height and time of field samples (e.g., Figs. 3-4) and elsewhere (Cooper et al., 1991). However, several other processes might contribute to generate variability, while others can act to dampen this variability. In the next section, we will discuss alternative causes, complementary and antagonistic, that contribute to the observed intraindividual δ xyl variances.

Fractionation at root or stem level
An increasing body of observations shows the occurrence of isotopic fractionation at the root level governed by root membrane transport (Lin and Sternberg, 1993;Vargas et al., 2017) or by unknown reasons (Zhao et al., 2016). Brinkmann et al. (2019) hypothesize that root level fractionation causes disparity when average RWU depth calculations based on δ 2 H X measurements are compared with those of δ 18 O X . However, it is difficult to imagine a scenario where root fractionation by itself can explain the observed diurnal fluctuations in δ xyl with height and time. Even if root fractionation significantly contributed to variation in δ xyl , we would still need to take into account diurnal fluctuation in RWU to explain the observed patterns. Isotopic enrichment of xylem water along the stem length was observed in association with stem transpiration (Barnard et al., 2006;Dawson and Ehleringer, 1993). However, this phenomenon is generally restricted to non-suberized plants and in woody branches in close vicinity to the evaporative surface of the plant (Dawson and Ehleringer, 1993). Isotopic enrichment can, therefore, not explain the variances in δ xyl observed in our empirical data, which were sampled within the main stem (data French Guiana) or from lignified branch segments distant from evaporative surfaces (data China and Germany).

Temporal and spatial soil dynamics
Soil water content can be extremely heterogeneous in the three spatial dimensions, as well as in time with complex dynamics of soil water movement. For example, hydraulic lift vertically redistributes soil water through the roots (Dawson and Ehleringer, 1993), which may change the water isotopic composition of the water mixture in the rhizosphere that is taken up by roots. Specifically, hydraulic lift redistributes and mixes the isotopic signal of depleted soil water in deeper layers with the enriched soil water signal in the rhizosphere in shallower layers. This should lead to lower variation in the soil water accessible to the plant, and hence less variation along plant height. Horizontal heterogeneity of water content may also affect δ xyl variance as soil water potentials and the isotope composition of soil water are interlinked. Under these conditions, it is important to understand how much the radial distribution of roots will naturally average out soil heterogeneity. However, note that heterogeneity in the soil does not automatically translate into variability in the xylem. Differential root water uptake driven by the diurnal fluctuation in water potential gradients in the soil-plant interface is still required to generate variability in the xylem isotopic signature.

Storage tissue and phloem enrichment
Storage tissues release water and sugars into the xylem conduits on a daily basis to support water transpiration demand (Goldstein et al., 1998;Morris et al., 2016;Secchi et al., 2017) or to repair embolism (Salleo et al., 2009;Secchi et al., 2017). Both water and sugars are transported in and out of storage tissue via symplastic pathways using plasmodesmata and aquaporins (Knipfer et al., 2016;Secchi et al., 2017), a pathway that has been linked to isotopic fractionation in roots (Ellsworth and Williams, 2007). Moreover, phloem transports photosynthetic assimilates that were produced in the leaves and are therefore potentially affected by transpiration fractionation (Gessler et al., 2013). Hence, these metabolic molecules might show higher values of δ 2 H and δ 18 O compared to RWU. Water release from storage or phloem tissue might locally alter δ xyl (White et al., 1985). Additionally, the time between water storage and release could bridge multiple days, and corresponding isotopic composition may reflect soil conditions preceding a dry spell when the isotopic signature of soil was less vertically stratified. It is evident that such dynamics are complex, and it is hard to predict how storage tissue and phloem enrichment affect observed δ xyl patterns. Importantly, xylem isotopic sampling cannot differentiate between water resulting from RWU or storage, and therefore we cannot exclude the possibility that tissue and phloem enrichment play a role. At a minimum this adds further uncertainty to RWU assessment. Water derived from storage tissues might also be present in larger fractions in higher parts of the plants, especially branches, as contamination accumulates as water moves upwards.
Unfortunately, to our best knowledge, empirical data on the isotopic composition of storage tissue and its spatiotemporal dynamics are absent in the literature. Future research should target the impact assessment of storage water on intraindividual δ xyl , allowing proper implementation in the model.

Diffusion processes
Diffusion is a process of net movement of molecules from a region of higher concentration to a region of lower concentration. Consequently, diffusion dampens δ xyl variability both in time and space within water xylem. Although the mutual diffusion coefficient of heavy water in normal water is very small and flow within vessels is laminar, other processes might still contribute to generating diffusion along the xylem. For example, as the water moves through the complex network of vessels, differences in velocities between vessels of different sizes cause some particles to move faster or slower than the average flow. According to the Hagen-Poiseuille law, the flow in each vessel is proportional to the fourth power of the vessel radius and the mean velocity is proportional to the square of the radius, thus potentially generating large differences in particle velocities depending on the vessel size distribution and other anatomical properties. Even within a single vessel, velocity is parabolic with a maximum flow velocity in the center and zero at the vessel walls.

A way forward
The observed large δ xyl variance and temporal dynamics in the empirical data suggest the need for a critical assessment of the stable isotope tracer technique for RWU studies. However, it also creates new opportunities. Since δ xyl variance and temporal dynamics herein likely relate to various plant physiological processes, the monitoring of variation in δ xyl can allow a more integrated understanding of plant water transport and hydraulic properties.
Combining a plant hydraulic model with in situ SF V , δ 2 H S,i and S,i,t can also help improve the robustness of RWU assessment and interpretation. Measurements of δ 2 H S,i and S,i,t at multiple depths, i.e., by installing soil water suction cups working in a vacuum (i.e., Rennenberg et al., 1996) and multiple soil matric potential sensors that measure at a high temporal frequency, should be especially valuable since the SWIFT model showed high sensitivity to alterations of this variable and these can be directly supplied as model inputs. At the same time, the availability of SF t measurements allows for the identification of the moment when water uptake from all root layers is at its maximum, which can be used to determine the optimal timing of sampling at a given height providing a more robust estimation of average RWU depth and uptake.
Alongside the modeling approach presented here, new ways to study δ 2 H X at a high temporal scale are strongly encouraged, for example, the pioneering work of Volkmann et al. (2016) in the development of an in situ continuous isotope measurement technique that offers the possibility for monitoring δ xyl at a sub-hourly resolution. This technique holds strong promise for further elucidating the natural δ 2 H X variances found within plants and the physiology processes from which these variances result. Such a high temporal resolution of isotope measurements, coupled with in situ monitoring of various environmental and plant biophysical metrics, is needed for both model improvement and further validation. Moreover, these seem inevitable to eventually differentiate all causal mechanisms of the observed intraindividual δ xyl variance.

Conclusions
A collection of empirical field data show pronounced variance and high temporal fluctuations in δ xyl . Moreover, these high temporal fluctuations in δ xyl emanate from basic plant hydraulic functioning, as model explorations show. We expect the observed δ xyl variance, and sub-daily fluctuations result, for a large part, from the mechanisms considered here, though various other physiological processes could also affect δ xyl .
Our theoretical explorations warn that variability in the isotope composition of plant xylem water can result in erroneous average RWU depth estimation and will complicate the interpretation and comparison of data; samples taken at different heights and times or plants differing in SF V may incorrectly show differences in average RWU depth. We further predict that various soil parameters and plant hydraulic parameters affect (i) the absolute size of the error and (ii) the probability of measuring δ xyl values that do not represent the well-mixed values during the plants' peak RWU. Hydraulic models, such as SWIFT, could help us to design more robust sampling regimes that enable improved comparisons between studied plants. We advocate the addition of SF t , which indirectly reflects diurnal RWU fluctuations, and S,i,t monitoring as a minimum in future RWU assessments since these parameters were predicted to be the predominant factors introducing variance in δ xyl from the SWIFT model exploration. However, soil texture and root permeability are also key measurements especially when comparing across species and sites.
Our findings do not exclude additional factors that impact the observed intraindividual δ xyl variance and temporal fluctuation as many processes can act simultaneously and are not mutually exclusive. Therefore, we strongly emphasize the need for more research. Directed studies that validate and quantify the relative impact of other plant physiological processes towards variance in δ xyl are a prerequisite before improved modeling tools can be developed. Shape parameter for the soil hydraulic properties (Clapp and Hornberger, 1978) Dimensionless B i The overall root length density distribution per unit of soil (not necessarily limited to the focal plant) m m −3 δ 2 H X,0,t Isotope composition of plant xylem water at stem base at time t In ‰ VSMOW δ 2 H X,h,t Isotope composition of plant xylem water at height h and time t In ‰ VSMOW δ 2 H S,i Isotope composition of soil water of the ith soil layer (constant over time) In ‰ VSMOW δ sample Isotope composition of water within a sample In ‰ VSMOW ˆ i,t Estimated water potential gradient between stem base and ith soil layer at time t derived from Eq. (8) m i,t Soil matric potential gradient between soil and roots at the ith soil layer at time t m H 2 O β 2 H X ; β 18 O X Normalized isotope composition of plant xylem water In ‰ VSMOW f i,t The fraction of water taken up in the ith soil layer at time t Dimensionless h Measurement height m i Soil layer index Dimensionless δ xyl Isotope composition of plant xylem water In ‰ VSMOW k i Soil-root conductance of the ith soil layer s −1 K max Maximum soil hydraulic conductivity m s −1 k R Effective root radial conductivity s −1 k S The conductance associated with the radial water flow between the soil and the root surface s −1 K S,i Soil hydraulic conductivity at the ith soil layer m s −1 l The approximated radial pathway length of water flow between bulk soil and root surface m LF Lumen fraction per unit sapwood area m 2 m −2 n Number of unique contributing water sources # sat Soil matric potential at soil saturation Delay before the isotope composition of xylem water at stem base reaches stem height h s θ sat Soil moisture content at soil saturation m 3 m −3 θ S,i,t Soil moisture content of the ith soil layer at time t m 3 m −3 z i Soil depth of the ith soil layer m Data availability. Both the French Guiana data and the SWIFT model are available on the GitHub repository HannesDeDeurwaerder/SWIFT. For the availability of the data collected in China and Germany, readers are referred to Zhao et al. (2014) and Magh et al. (2020), respectively.
Author contributions. HV, MDV, and PB supervised and provided guidance throughout all aspects of the research. HDD, MDV, and HV designed the study. HDD, KK, RKM, JDM, LW, and LZ collected and processed the empirical datasets. The model was developed and coded by HDD, MDV, MD, and FM. All authors contributed to the interpretation of the results and the text of the paper.