Isotopic approaches to quantifying root water uptake and redistribution: a review and comparison of methods

. Plant root water uptake (RWU) and release (i.e., hydraulic redistribution – HR, and its particular case hydraulic lift – HL) have been documented for the past five decades from water stable isotopic analysis. By comparing the (hydrogen 10 or oxygen) stable isotopic composition of plant xylem water to those of potential contributive water sources (e.g., water from different soil layers, groundwater, water from recent precipitation or from a nearby stream) authors could determine the relative contributions of these water sources to RWU. Other authors have confirmed the existence of HR and HL from the isotopic analysis of the plant xylem water following a labelling pulse. In this paper, the different methods used for locating / quantifying relative contributions of water sources to RWU (i.e., 15 graphical inference, statistical (e.g., Bayesian) multi-source linear mixing models) are reviewed with emphasis on their respective advantages and drawbacks. The graphical and statistical methods are tested against a physically based analytical RWU model during a series of virtual experiments differing in the depth of the groundwater table, the soil surface water status, and the plant transpiration rate value. The benchmarking of these methods illustrates the limitations of the graphical and statistical methods (e.g., their inability to locate or quantify HR) while it underlines the performance of one Bayesian 20 mixing model, but only when the number of considered water sources in the soil is the highest to closely reflect the vertical distribution of the soil water isotopic composition. The simplest two end-member mixing model is also successfully tested when all possible sources in the soil can be identified to define the two end-members and compute their isotopic compositions. Finally, future challenges in studying RWU with stable isotopic analysis are evocated with focus on new isotopic monitoring methods and sampling strategies, and on the implementation of isotope transport in physically based 25 RWU models.


INTRODUCTION
Root water uptake (RWU) is defined as the amount of water abstracted by a root system from soil over a certain period of time. Understanding the relation between the distributions of soil water, roots, RWU location and magnitude, and root hydraulic properties is important for managing soil water and plant water status (e.g., by irrigation), developing new plant genotypes more tolerant to drought or tackling ecological questions in water-limited ecosystems, such as the competition for 5 soil water by different plants.
RWU is principally driven by evaporative flux taking place in the leaves (i.e., transpiration) and its magnitude depends on the atmospheric evaporative demand and stomatal opening. The latter depends amongst others on leaf water status and stress hormonal signals from the roots transported to the leaves (e.g., Huber et al., 2015;Tardieu and Davies, 1993). Leaf water status and hormonal signals are related to the soil water potential distribution and to the plant hydraulic architecture (Huber 10 et al., 2015). The spatial distribution of RWU is very variable in time in space, depends on the presence of roots but also on the root's ability to extract water. This ability is a function of radial conductivity but axial conductance may also limit water flow in younger roots or when cavitation occurs. The flux of water depends also on soil water availability: a highly conductive root segment will not be able to extract water from a dry soil. Locally, this is the difference of water potential between the root and the soil which drives RWU, and its magnitude is controlled by the radial hydraulic resistances in the 15 rhizosphere, at the soil root interface and in the root system (Steudle and Peterson, 1998). The actual RWU profile is thus a combination of different aspects: the root ability to extract water (characterized by the amount of roots and their hydraulic properties), the ability of the soil to fulfill the plant water demand, and the water potential difference between soil and root (Couvreur et al., 2014).
Plants have numerous mechanisms to cope with heterogeneous soil water distribution, and adapt their RWU rate distribution: 20 adaptive root growth, adaptive root conductivity (Javaux et al., 2013), exudation (Carminati et al., 2016). A particular process, which has attracted the attention of plant breeders and ecologists is the ability of plants to extract water from non or less water limited soil areas with potentially low root lengths densities, known as root water uptake compensation (Heinen, 2014). To describe the RWU rate in soils, we will use the root water uptake flow per volume of soil, defined as S [L 3 L -3 T -1 ] by reference to the sink term of the Richardson (1922) equation . According to Couvreur et al. 25 (2012), root compensation is defined in the present article as the process that decreases or increases RWU at a certain location compared to the water uptake from that location when the soil water potential would be uniform in the root zone.
Thus, the distribution of the S(x,y,z) is a sum of two spatially distributed components: where x, y and z are the 3-D spatial coordinates, S uniH is a term proportional to the root distribution and S comp the 30 compensatory part of the RWU distribution. The first term on the right-hand side of Eq. (1) is always positive while the second one can be either positive or negative. Figure 1 illustrates how this equation affects S distribution in a onedimensional (1D) space. When there is no compensation (S comp (x,y,z) = 0), the RWU distribution follows the root distribution Biogeosciences Discuss., doi: 10.5194/bg-2016-410, 2016 Manuscript under review for journal Biogeosciences Published: 5 October 2016 c Author(s) 2016. CC-BY 3.0 License.
(i.e., highest at the surface and lowest in the deepest layer, Fig. 1a). When S comp (x,y,z) < 0 but its absolute value is lower than S uniH (x,y,z), then S(x,y,z) is positive and different from the root vertical distribution. In case S uniH (x,y,z) is small, as in Fig. 1c, S comp (x,y,z) can locally be higher in absolute value and S(x,y,z) can be locally negative which implies that there is a water efflux out of the root.
The water efflux at certain locations is called root hydraulic redistribution (HR, Burgess et al., 1998) or hydraulic lift (HL, 5 Richards and Caldwell, 1987) as a specific case of HR in which fluxes in the root system are vertically upward. In their review, Neumann and Cardon (2012) discussed that the magnitude of HR observed in different studies varied from 0.03 mm d -1 (brasilian Cerrado, Scholz et al., 2010) to 3.50 mm d -1 (Artemisia Tridenta, Ryel et al., 2003). Several authors have also raised the question of the "ecohydrological interest" for a plant to release water to the upper/dryer soil layers, therefore potentially providing water to shallow-rooted plants and enhancing competition for space and nutrients. Some studies 10 suggested that HL could increase nutrient mobility and enhance biogeochemical processes by providing moisture to the dryer soil layers (Caldwell et al., 1998;Snyder et al., 2008).
Despite its importance, there is a lack of measurements of RWU, related to the difficulty of measuring root and soil water fluxes. Often soil water content change is used as a proxy for RWU. Yet, as change of soil water content with time is not due to root extraction only (i.e., soil water redistribution can also occur), the assessment of RWU based on water content 15 distribution alone is not possible in conductive soils . Rather, the full soil water flow equation accounting for root uptake and soil water redistribution must be solved in an inverse mode, and, with an accurate knowledge of soil and root properties RWU distribution can be inferred (Guderle and Hildebrandt, 2015;Hupet et al., 2002;Musters and Bouten, 1999;Vandoorne et al., 2012). Nuclear magnetic resonance (NMR) imaging has been suggested as an adequate technique to measure water flow velocity in xylem vessels but no application exists yet on living roots in soils (Scheenen et 20 al., 2000). More recently, Zarebanadkouki et al. (2012) could measure for the first time RWU in porous media by combining a tracer experiment monitored by neutron tomography with inverse modelling of a transport equation. Yet, this was done under controlled conditions while there is no standard method to monitor three dimensional water uptake distribution of growing roots in situ. In woody plants, in which roots are thick enough, Nadezhdina et al. (2010; used sap flow measurements in roots to quantify hydraulic redistribution. 25 Since the seminal work of Zimmermann et al. (1967) which reported that RWU of Tradescantia fluminensis occurred in the absence of fractionation against water oxygen stable isotope, water stable isotopologues ( 1 H 2 H 16 O and 1 H 2 18 O) have been frequently used to identify and quantify root water uptake and redistribution in soils through the measurements of their natural (and artificial) isotopic abundances. Methods include simple graphical inference to more sophisticated statistical methods, i.e., two-end members and multi-source linear mixing models. While the former attempts to locate the "mean root 30 water uptake" in the soil, the latter category of methods provides profiles of relative contributions to transpiration flux across a number of defined soil layers.
This present paper (i) aims at reviewing these methods and (ii) proposes to compare them against each-other during a series of virtual experiments differing in the water and isotopic statuses in the soil and the plant. Prior to the review and intercomparison, the paper reports on the mechanisms at the origin of the spatiotemporal dynamics of natural isotopic abundances in soil and on the background knowledge of isotopic transfer of soil water to and from roots. Finally, we address future challenges to be undertaken such as the dynamic isotopic assessment of HR. We also evoke opportunities offered by novel 5 isotopic monitoring tools which provide unpreceded high frequency isotopic measurements.

THEORETICAL BACKGROUNDS
The temporal and spatial variations in natural isotopic abundances observed within the soil-plant system allow for reconstruction of S profiles. These variations result from isotope-specific fractionation between different phases at thermodynamic equilibrium and during non-equilibrium phase transition, i.e., when there is a net flux between different 10 phases as for instance during evaporation. In this section, we briefly review process based analytical models accounting for isotopic fractionation that were first proposed for (i) free water (section §2.1) and (ii) later for matric-bound water in a bare soil (section §2.2). Finally, (iii) we report on the absence of isotopic fractionation during RWU for most of the documented plant species and on the simple mixing model which is at the basis of any isotopic study on RWU (section §2.3).

Isotopic effects during free water evaporation 15
In a closed liquid water-water vapor isothermal system at water vapor saturation (relative humidity = 100 %) or at thermodynamic equilibrium, the difference between the liquid and vapor (hydrogen or oxygen) isotopic compositions (δ l and δ v [-, expressed in ‰ relative to the Vienna-Standard Median Ocean Water international isotope reference scale], Gonfiantini, 1978) is a function of the system temperature solely and is named "equilibrium isotopic fractionation factors" (α eq [-], for a complete list of symbols see Appendix A). Majoube (1971) and Horita and Wesolowski (1994), among other 20 authors, gave empirical expressions (i.e., closed-form temperature dependent equations) for these equilibrium fractionation factors.
When the system is no longer closed and a difference in water vapor partial pressure exists between the air layer in direct contact with the liquid surface and the atmosphere above (referred to as "free" atmosphere), water vapor is transferred from the liquid phase to the air layer, i.e., evaporation (E) occurs. In analogy to an electrical circuit (i.e., a Rideal-Langmuir 25 linear-resistance model, Brutsaert, 1982), E can be calculated from the vapor pressure difference and a transfer resistance (r) to vapor transport across the air layer between the evaporating surface and the free atmosphere. Following the same electrical analogy, vapor transport of isotopologues (E i ) is a function of the difference of vapor isotopic composition between the air layer in contact with liquid water and the free atmosphere and a transfer resistance (r i ). In a simple yet comprehensive model Craig and Gordon (1965) divided the air layer into two consecutive layers with different aerodynamic conditions. In a first 30 Biogeosciences Discuss., doi:10.5194/bg-2016-410, 2016 Manuscript under review for journal Biogeosciences Published: 5 October 2016 c Author(s) 2016. CC-BY 3.0 License.
sub-layer (with transfer resistances r diff and r diff,i for lighter and heavier isotopologues, respectively), vapor transfer is purely diffusive whereas in a second layer (with transfer resistances r ad and r ad,i for lighter and heavier isotopologues, respectively) it is purely advective, i.e., controlled by turbulence. Since the diffusion coefficient of the heavier isotopologues are smaller than that of 1 H 2 16 O, r i is greater than r. Note that the difference between r i and r originates from the difference between r diff and r diff,i whereas r ad = r ad,i . In this model, thermodynamic equilibrium conditions still prevail in the air layer in contact with 5 the evaporating surface, termed "liquid-vapor interface" of isotopic composition δ l-v . The isotopic composition of the evaporated water vapor (δ E ), defined as the ratio E i /E, depends on the ratio of the resistances r i /r. The latter ratio, named "kinetic isotopic fractionation factor" (α K , [-]) depends on the relative importance (or development) of each sub-layer and contributes producing an evaporated vapor depleted in heavy isotopologues with respect to the vapor at the liquid-vapor interface (i.e., δ E < δ l-v ). In turn, and depending on the turnover of the system (ratio E/V with V being the volume of the 10 evaporating liquid), but also on the evaporation state (i.e., permanent or transient), the liquid phase enriches itself in the heavy isotopologues. Finally, when both E and V are constant over time, meaning that the loss of water is compensated by a source of constant isotopic composition δ source , an "isotopic steady state" might be reached where, by mass balance, δ E = δ source . For a thorough review of the evaporation model of Craig and Gordon (1965), the reader is referred to Gat (1996) and to the more recent paper by Horita et al. (2008). 15 In a two-dimensional (δ 18 O, δ 2 H) space, meteoric waters (e.g., precipitation, river water, groundwater) formed by equilibrium processes (i.e., condensation of water vapor) fall onto a line whose slope equals approximately eight and whose theoretical value is the ratio at the temperature of condensation. On the other hand, the water vapor produced during a non-equilibrium process, such as evaporation, fall onto a so-called "evaporation line" with a slope of generally lower than six and greater than 2. This is explained by the fact that  Gat (1971) showed that the value of this slope was fairly approximated by the Craig and Gordon (1965) model, which was recently tested by Rothfuss et al. (2015).

Isotopic effects during bare soil evaporation and leaf transpiration
The Craig and Gordon (1965) model, originally developed for free evaporating water was later adapted to bound-to-matrix soil water. In a study that laid the basis for future work in isotopic ecohydrology, Zimmerman et al (1967) provided a steady-25 state analytical solution for soil water isotopic composition (δ S ) in a water-saturated isothermal bare sand profile from which water evaporated at a constant rate. Under these steady-state and isothermal conditions, the upward (convective) liquid flux of isotopologues, triggered by evaporation and rising from deeper layers equals the downward (diffusive) isotopic flux from the evaporating surface which is enriched in the heavy stable isotopologues due to evaporation. A profile is obtained (Fig.   2a, black line) whose exponential shape depends on boundary conditions, i.e., the source water (e.g., groundwater) and 30 surface water isotopic compositions (δ source and δ surf ), the diffusion coefficient of the isotopologues in water, and of a soil "tortuosity factor", conceptually defined as the ratio of the geometrical to actual water transport distance. Allison Biogeosciences Discuss., doi:10.5194/bg-2016-410, 2016 Manuscript under review for journal Biogeosciences Published: 5 October 2016 c Author(s) 2016. CC-BY 3.0 License.
(1983) extended this formulation to a non-saturated sand column evaporating at isotopic steady state (δ E = δ source ). In this case, the evaporating surface (i.e., the liquid-vapor interface) can be located below the soil surface and splits the profile into two regions where isotopic transport predominantly occurs either in the vapor phase above or in the liquid phase below it. At isotopic steady state, the maximal isotopic enrichment is at the evaporation front (δ EF at soil depth z EF ) and can be simulated with the Craig and Gordon (1965) model. The isotopic composition of the soil residual adsorbed water in the "vapor region" 5 above the evaporation front can be obtained by assuming thermodynamic equilibrium conditions and by applying Fick's law, and is shown to decrease linearly towards the value of the liquid water at the soil surface which is at thermodynamic equilibrium with the ambient atmospheric water vapor (Fig. 2a, gray line). Finally, note that at transient state (δ E ≠ δ source ), the maximal isotopic enrichment in the soil profile does not point to the location of the evaporation front as was demonstrated by Rothfuss et al. (2015). Instead, the depth where the steepest gradient in the isotopic profile is observed 10 corresponds to the evaporation front.
In a two-dimensional (δ 18 O, δ 2 H) space, liquid soil water sampled below the evaporation front will plot on an "evaporation line" with a slope typically lower than six and greater than two, depending on atmospheric and isotopic forcing, as a result of kinetic processes during evaporation. Above the evaporation front and at isotopic steady-state, soil liquid water is in equilibrium with a mixture of atmospheric water vapor (δ 18 O-δ 2 H slope ~8) and evaporated soil water vapor rising from the 15 evaporation front (2 < δ 18 O-δ 2 H slope < 6) ( Barnes and Allison, 1988;Brunel et al., 1995;DePaolo et al., 2004). As a result, an intermediate value for the slope is expected, depending on the mixing ratio of atmospheric water vapor to evaporated soil vapor at a given soil depth.

Isotopic transfer to and from roots
As opposed to the removal of water vapor by evaporation, RWU has been described in a number of studies and over a wide 20 variety of plant species not to be associated with (kinetic) isotopic fractionation (Bariac et al., 1994;Dawson and Ehleringer, 1993;Thorburn et al., 1993;Walker and Richardson, 1991;Washburn and Smith, 1934;White et al., 1985;Zimmermann et al., 1967). Consequently, for plants growing in homogeneous external conditions, e.g., in hydroponic solution, root xylem sap water and external water have the same isotopic compositions. In natural soils where the liquid phase is not homogeneous and a vertical gradient of isotopic composition due to evaporation exists, the root system takes up water at 25 different depths having thus different isotopic compositions.
Assuming that water transport time in roots is negligible, the isotopic concentration of the xylem sap water at the root tiller with δ S (x,y,z) [-, expressed in ‰] the isotopic compositions of soil water at coordinates (x,y,z). Mostly, a one dimensional description of root water uptake is used assuming that  S and RWU do not vary in the horizontal direction and δ S is obtained for discrete soil layers of depths z j (j ∈ [1,n]) and thickness Δz j = z j+1 -z j . It is usually further hypothesized that J Ti equals the transpiration flux T [L 3 T -1 ] (low to no plant capacitance or phloem-xylem contact): 10 where q Ti = J Ti /(x  y) = T/(x  y) represents the sap flow rate in the root tiller per unit surface area [L T -1 ].
δ Ti can be accessed at different locations in the plant depending on the species, but the sampling location should not be affected by evaporative enrichment in heavier isotopologues or back-diffusion of the isotopic excess accumulated at the sites of transpiration (stomatal chambers) in the leaf. For grasses and nonwoody plants, this is done by sampling the root crown 15 (e.g., Leroux et al., 1995), the aerial nodal roots (e.g., Asbjornsen et al., 2007), the meristematic petiole, or else the collars (e.g., tillers) at the base of the plant (e.g., Dawson and Pate, 1996;Sánchez-Perez et al., 2008). In the case of ligneous plants the fully suberized stem (Asbjornsen et al., 2007) or sapwood (e.g., White et al., 1985) is sampled. On the other hand, δ S is usually measured by sampling soil profiles destructively. Finally, water from plant and soil is predominantly extracted by cryogenic vacuum distillation (Araguás-Araguás et al., 1995;Ingraham and Shadel, 1992;Koeniger et al., 2011;Orlowski et 20 al., 2013;West et al., 2006). Lin and Sternberg (1992) and Ellsworth and Williams (2007), amongst other authors, reported however that for some xerophyte (plants adapted to arid environments, e.g., Prosopis velutina Woot.) and halophytes species (plants adapted to saline environments, e.g., Conocarpus erecta L.), and mangrove species (e.g., Laguncularia racemosa Gaert.), RWU led to fractionation of water hydrogen isotopologues. For mangrove species, it was hypothesized that the highly developed 25 Casparian strip of the root endodermis would force water moving symplastically (i.e., inside the cells) and therefore crossing cell membranes (Ellsworth and Williams, 2007 tranport and leads to a situation where xylem sap water is depleted in this isotopologue with respect to source water.
Meanwhile, this affects to a much lesser extent 1 H 2

18
O transport, so that no detectable isotopic fractionation of water oxygen isotopologues is observed. It can be concluded that, for the majority of the studied plant species, either RWU does not lead to isotopic fractionation or its magnitude is too low to be observable.
Hydraulic redistribution (e.g., hydraulic lift) can be conceptualized as a reverse RWU, defined as a negative S. In such case, 5 Eq. (4b) should only account for the positive S. It can be done by assuming that in this equation δ S is 0 when S < 0.
Finally, plant water samples will, similarly to soil water samples, also fall onto an "evaporation line" of a slope lower than eight in a two-dimensional (δ 18 O, δ 2 H) space .

METHODS FOR CHARACTERIZING RWU FROM STABLE ISOTOPIC ANALYSIS
We distinguish two classes of methods: (i) the graphical method for inferring the "mean root water uptake depth" ( z [L]) 10 ( §3.1), and (ii) statistical methods based on end-member mixing analysis (EMMA) (Barthold et al., 2011;Christophersen and Hooper, 1992) for identifying x j [-], the contribution to RWU of some plant water source j (e.g., water in some soil layer, groundwater, water from recent precipitation, or else from a nearby stream, §3.2). All methods have in common to use an inverse modeling approach: the RWU distribution is obtained by optimizing model input parameters until the simulated δ Ti and/or the simulated soil isotopic profiles fit to the isotopic measurements. 15 Table 1 summarizes the 21 isotopic studies reviewed in this paper that use either one of the two classes of methods.

Graphical inference (GI)
This straightforward approach defines the "mean root water uptake depth" z , as the depth where  S = δ Ti . z conceptually indicates to the soil depth where the plant root system, represented as one unique root, would extract water from.
There are cases where z cannot be unambiguously identified (e.g., 1 z and 2 z of case 2, Fig. 2b) due to the non-monotonic 20 character of the δ S profile (shown in black dashed line, case 2 of Fig. 2b). In order to define a mean RWU depth for such a case one can derive a monotonously decreasing δ S profile by smoothing the profile (shown as symbols in Fig. 2b), e.g., by averaging  S in a number of layers using the following mass balance: where J represents the set of depths that belong to the J th soil layer, with θ [L 3 L -3 ] and Δz j [L] the soil volumetric water 25 content and thickness of the soil layer centered around depth z j . Due to this smoothing, the vertical resolution may be drastically reduced. In the example presented Fig. 2b where a uniform θ profile is assumed, the δ S,J profile intersects with the vertical line of value δ Ti deeper than for the initially non-monotonic δ S profile, i.e., z (case 2, integrated δ S profile) < . Some authors rule out solutions in case of multiple mean root water uptake depths, e.g., by excluding the z solutions where soil water content was low and/or soil water potential was high in absolute value (e.g., Li et al., 2007; see Table 1).
Note that while Eq. (5) provides a representative value for the isotopic composition that would be measured in soil layer J as a function of those of the water in the set of depths, δ S,J is however equivalent to the isotopic composition "sensed by the plant" only if the root profile is homogeneous, i.e., when RLD is constant with depth in that particular soil layer J. 5 The method of graphical inference may not only provide z but also its uncertainty caused by the uncertainty in measuring δ Ti (e.g., based on the precision of the isotopic analysis and/or sampling natural variability, shown as gray stripe in Fig. 2b). The steeper the soil water isotopic profile, the larger is the uncertainty in determining z is. Figure 2b illustrates this with estimated minimum and maximum z for the monotonic δ S profile and for the vertically averaged profile. In the latter case, the possible range of z is the largest. These ranges give first quantitative indication of variance around z . Finally, for a 10 complete "graphical assessment" of the variance of z , one should also consider the uncertainty associated with measurements of the δ S profile (not shown).

Two end-member (TM) mixing model
The TM method is a particular case of end-member mixing analysis (EMMA) and is based on the concept that (i) a plant 15 extracts water from two predominant water sources A and B (e.g., water in distinct upper and lower soil layers, or groundwater and recent precipitation water etc.) in given proportions, (ii) there is no isotopic fractionation during water uptake, and (iii) there is a complete mixing inside the plant of the contributing water sources A and B to RWU. The mass conservation for isotopologues gives: In this approach, δ Ti is therefore defined as the mean value of the isotopic compositions of water sources A and B (δ A and δ B ) weighted by the proportions to J Ti of water volume extracted by the plant from water sources A and B, i.e.,  i.e., the analytical precision of the isotopic analyser (e.g., isotope ratio mass spectrometer, laser-based spectrometer), or by 5 taking into account additional errors involved with sampling procedure and vacuum distillation technique (see e.g., Rothfuss et al., 2010). Equation (8b) also shows that, independently of the values considered for , indicating that the two end-members should have as much as possible distinct isotopic compositions for a low standard error of x. Therefore, it is especially important, e.g., for partitioning between water from an upper and lower portion of the soil profile, to properly define the thickness of these layers, so that they have distinct isotopic 10 compositions, and that the difference is considerably larger than the precision of the isotopic measurements. Figure 3 shows for example that when (i) x is evaluated at 10 % and (ii) are estimated being equal to 0.02 ‰ (dark blue solid line), (δ A -δ B ) should be greater than 0.75 ‰ (in absolute term) in order to reach a σ x value lower than 5 %, i.e., more than 37 times the error made on δ A , δ B , and δ Ti . To obtain the same standard error for x in case of a higher standard error on the estimation of δ A , δ B , and δ Ti ) ‰ 1 . 0 and , , (e.g.,  Figure 3). This certainly highlights the advantage of artificially labelling soil water with water enriched (or depleted) in heavy isotopologues for a more precise assessment of the relative contribution of soil water sources to RWU, as mentioned by Moreira et al. (2000). In another study, Bachmann et al. (2015) labeled the upper and lower portion of the soil profile in a natural temperate grassland with 18 O-enriched and 2 H-enriched water, respectively. They defined two distinct (upper and lower) soil water 20 sources, for which they calculated the corresponding δ 2 H or δ 18 O on the basis of measured soil water isotopic profiles and using Eq. (5). They could find evidence against the so-called hypothesis of "niche complementarity" regarding plant water use, which states that RWU of competitive plant species is spatially and temporally distinct, and that this distinction is stronger at high species richness. Figure 3 illustrates also that for given (δ A -δ B ), value" for a low σ x is 50% (showed by the orange lines). 25 Table 1 displays a sample of studies that used the two end-member mixing approach. Authors could distinguish between uptake of irrigation and precipitation water (Goebel et al., 2015), precipitation and groundwater (White et al., 1985), soil water and groundwater (McCole and Stern, 2007), or else between stream water and soil water (Dawson and Ehleringer, 1991;McDonnell, 2014). Thorburn and Ehleringer (1995) could for instance locate the dominant source for RWU, i.e., Biogeosciences Discuss., doi: 10.5194/bg-2016-410, 2016 Manuscript under review for journal Biogeosciences Published: 5 October 2016 c Author(s) 2016. CC-BY 3.0 License. groundwater for their mountain and floodplain test-site and water from the soil between 0.3 and 0.4 meters depth for their cold desert test-site. Other authors (e.g., Brunel et al., 1995) combined two mixing equations, i.e., one for each isotopologue, into a single one. As infrared laser-based spectrometry now enables simultaneous measurements of δ 18 O and δ 2 H at lower cost, we believe that this dual-isotope approach (referred as "D" in Table 1) will or should gain in importance in isotopic studies, especially in the context of pulse labelling experiments, which can "disconnect" the strong correlation between soil 5 water δ 18 O and δ 2 H, therefore provide two independent mixing equations, one for each isotopologue.

Multi-source (MS) mixing models
When there are more than two identified plant water sources contributing to RWU, e.g., water from different layers j (j ∈ [1, N]) in soil the profile, Eq. (7) becomes: with N the number of plant water sources (e.g., soil layers) and    N j j x 1 1 . As there are more water sources than (number of mixing equations + 1), there is not a unique solution but an infinite range of possible solutions. However, some of these solutions are not likely or possible based on background information or knowledge. A range of solutions that is most likely based on prior information can be obtained using Bayesian methods. In the method proposed by Phillips and Gregg (2003), the isotopic composition calculated for each considered x j combination (δ Ti ) is compared with the measured value (δ Ti,m ). The 15 number of combinations depends on the value of contribution increment (i, %, typically 5 or 10 %) and the combinations for which δ Ti meets the following requirement are selected: where τ [-, expressed in ‰], standing for "tolerance", usually accounts for precision of the isotopic measurements or possible errors during sampling and vacuum distillation steps. This multi-source mixing model approach strongly depends on τ and i, 20 which therefore should be carefully chosen by the user. A smaller i also refines the analysis. For this, the program "IsoSource" (https://www.epa.gov/sites/production/files/2015-11/isosourcev1_3_1.zip) is available (Phillips et al., 2005). Wang et al. (2010) compared the outcome of the GI and MS approaches and came to the conclusion that even though the latter did not solve the non-uniqueness problem and provided diffuse patterns of frequency that were difficult to interpret in some cases (e.g., in case of a non-monotonic isotopic profile), it had the advantage over the former method of providing a 25 systematic and quantitative assessment of ranges of relative contributions. Parnell et al. (2010) proposed to overcome two limitations of the approach of Phillips and Gregg (2003), i.e., its inability to (i) account for uncertainty in the estimations of δ Ti and of the water sources isotopic compositions δ S,j , and (ii) provide a optimal solution rather than ranges of feasible solutions. For this they use a Bayesian framework (for details see also Erhardt Biogeosciences Discuss., doi:10.5194/bg-2016-410, 2016 Manuscript under review for journal Biogeosciences Published: 5 October 2016 c Author(s) 2016. CC-BY 3.0 License. and Bedrick, 2013; Moore and Semmens, 2008;Parnell et al., 2013), which allows uncertainty in the x j proportions and incorporates a residual error term ε j (normally distributed with mean equal to zero and variance σ 2 ): Note that the terms of (i) trophic enrichment factor (TEF [-, expressed in ‰], see, e.g., meta-analysis of Vanderklift and Ponsard, 2003) and (ii) isotope concentration dependency Phillips and Koch, 2002) originally 5 incorporated in the formulation of Parnell et al. (2010) for other applications are not present in Eq. (9') since (i) no isotopic fractionation during RWU is assumed and (ii) isotope concentration dependency applies only for situations where isotopic compositions of different elements are measured and available. Parnell et al. (2010) developed the program "Stable Isotope Analysis in R" (SIAR, https://cran.rproject.org/src/contrib/siar_4.2.tar.gz) in which the initial (a priori) x j distribution is by default the Dirichlet distribution, of 10 which information can be partly specified by the user. A posteriori x j distribution is obtained by fitting the linear model to data via a Metropolis-Hasting (Hastings, 1970;Metropolis et al., 1953) Markov Chain Monte Carlo algorithm. Prechsl et al. (2015) apply both graphical and Bayesian approaches to evaluate the shift in z and change of RWU profile following drought treatments (approx. 20 to 40 % precipitation reduction with transparent rainout shelters) in both extensively and intensively managed grasslands. From both approaches it appeared that a shift in z was inexistent or not 15 observable from isotopic analyses. Another recent application of the Bayesian approach was performed by Volkmann et al.
(2016b), who took advantage of a newly developed soil isotopic monitoring method to confront high frequency δ S profiles time series to time series of δ Ti (indirectly obtained from the isotopic measurement of the transpired water vapor and assuming isotopic steady state, i.e., δ Ti = δ T ) following a labelling pulse (see Table 1 for details on the study).

INTER-COMPARISON OF METHODS 20
We tested and compared the different methods (GI, TM, MS) during a series of virtual experiments. Mean RWU depths (provided by the GI method) and x j distribution (provided by the TM and MS methods) were determined from soil and xylem water oxygen isotopic composition distributions. While the former information was prescribed to the different methods, the latter was calculated with the physically based analytical RWU model (referred to as AM) of Couvreur et al. (2012) (see Appendix B1 for a description of the model and Appendix B2 on how it was run for the inter-comparison). 25

Scenario definition
We developed eight virtual plausible scenarios of soil-plant systems under different environmental conditions. Each environmental condition was defined as a combination of different total soil water potential distributions (resulting from the Biogeosciences Discuss., doi: 10.5194/bg-2016-410, 2016 Manuscript under review for journal Biogeosciences Published: 5 October 2016 c Author(s) 2016. CC-BY 3.0 License. location of the groundwater table and weather conditions), soil water oxygen isotopic composition profiles, and actual transpiration rate. The groundwater table was either shallow at -1.25 m depth (prefix "Sh") or deep at -6 m depth (prefix "De"); the soil water potential was considered to be at static equilibrium below the groundwater level; the soil surface was either dry under evaporative conditions (suffix "Dr"), or wet, e.g., shortly after a rain event (suffix "We"); the transpiration rate was either low (e.g., relevant at night, T = 0.01 mm h -1 , suffix "_lT") or high (T = 0.30 mm h -1 , suffix "_hT"). They all 5 relied on a common measured root length density vertical distribution of Festuca arundinacea. Table 2 reports the input data. Note that, as hypothetized in Eq. (4b), transpiration and sap flow rates (i.e., per unit of surface area [L T -1 ]) were considered as equal.

Setup of the models
The two end-member mixing approach (TM) was tested against the isotopic data for two different cases: (i) two conjoint soil 10 layers spreading from 0 -0.225 m and 0.225 -2.00 m and (ii) two disjoint soil layers spreading from 0 -0.225 m and 1.75 -2.00 m. The latter case was designed to evaluate the impact of lacunar soil isotopic information on the calculation of x, i.e., when not all potential water sources are properly identified. Representative values of water oxygen isotopic compositions for these soil layers (δ S,J , J ∈ [I,II]) were obtained from the mass balance (Eq. (5)) after interpolation of the measured soil water content and δ S profiles at a 0.01 m vertical resolution. 15 For the multi-source mixing approaches of Phillips and Gregg (2003) (MSPG) and Parnell et al. (2013) (MSPa), the number of potential water sources was initially fixed to three, i.e., water from the soil layers I (0.000-0.050 m), II (0.050-0.225 m), and III (0.225-2.000 m). Upper and lower boundaries of these layers were defined to reflect the exponentially shaped (monotonic) δ S profiles (experiments ShDr and DeDr) or to smooth the non-monotonic δ S profiles observed during experiments ShWe and DeWe. MSPG and MSPa were also tested for eight soil layers (i.e., as many layers as measurement 20 points, I: 0.000-0.020, II: 0.020-0.050, III: 0.050-0.110, IV: 0.110-0.225, V: 0.225-0.400; VI: 0.400-0.750, VII: 0.750-1.500, and VIII: 1.500-2.000 m). Increment and tolerance of the MSPG method were fixed at 10 % and 0.25 ‰, respectively.
Similarly to the TM approach, profiles of δ S,J (J ∈ [I,III] or [I,VIII]) were obtained from the mass balance (Eq. (5)) after interpolation of the measured soil water content and δ S profiles at a 0.01 m vertical resolution.
Finally for the MSPa method, uncertainty associated with δ S measurements was set to 0.2 ‰ and the function 25 siarmcmcdirichletv4 of the SIAR R-package was run 500000 times (of which 50000 runs where discarded).
For a detailed description of the inter-comparison methodology, refer to Appendix C. description on how uncertainty was assessed, refer to Appendix C). In general, at high T the compensation was negligible and the S AM profile was mainly proportional to the RLD profile (Fig. 4b, d, f, and h). The only exception was a soil with deep groundwater table and dry surface, where this dry layer limited root water uptake (DeDr_hT). At lower transpiration demand, the S profile predicted by the Couvreur et al. (2012) model generally differed from the RLD profile (Fig.4a, c , e, and g) due to the fact that the second term of Eq. (1) (i.e., S comp , see also Eq. (B4) and (B4') in Appendix B) was 5 proportionally larger. Water uptake from the upper layer was always more than proportional to the RLD, when this layer was wetter, and vice versa. Water release to the soil (i.e., HR) was observed only for the soil with the deep groundwater table and dry upper layer (DeDr_hT, Fig. 4e). From the graphical method GI, either a single or two distinct solutions for z (displayed as gray-shaded horizontal stripes) could be retrieved, depending on the monotonic/non-monotonic character of the δ S profile, and ranged between -0.02 and -0.95 m. 10  Figure 5) were due to the small difference between the isotopic compositions of the defined soil water sources δ S,I (-6.0 ‰) and δ S,II (-5.3 ‰) as illustrated in section §3.2.1. Figure 5b gives the relative contribution to T of the layer 1.75 -2.00 m in case of two disjoint soil layers, i.e., when not all potential water sources are accounted for into the calculation of δ S,I and 20 δ S,II . In this case there were important disparities between x TM and x AM . The mean absolute difference between these two estimates was equal to 43,5 (±17.8) %. Omitting some of the potential water sources contributing to T had in this second case the consequence of artificially overestimate the contribution of the lowest layer. We therefore suggest to always attempting to fully characterize the soil isotopic profile before aggregating the isotopic information when defining the two water sources. 25 Figure 6 gives the relative contributions from soil layers I, II, and III (upper, middle, and lower panel, respectively) to T following the method of Phillips and Gregg (2003)  x III_AM refer to their uncertainty by accounting for the uncertainty of the input data. As for Fig. 5, δ Ti_AM is reported above 30 each plot along with its standard deviation. x J_MSPG probability distribution was observed to be either narrow (e.g., DeDr_lT / layer I, Fig. 5m) or broad (e.g., DeWe_hT / layer I), i.e., the range of the possible solutions for x J_MSPG was relatively small Biogeosciences Discuss., doi: 10.5194/bg-2016-410, 2016 Manuscript under review for journal Biogeosciences Published: 5 October 2016 c Author(s) 2016. CC-BY 3.0 License. or large (10 and 100 % respectively for these two examples). In general, both MSPG and MSPa statistical methods agreed well with each other: the x MSPa most frequent value (MFV, at the peak of the density distribution curve) was in most case either located near the median value of the x MSPG probability range (e.g., ShWe_lT / layer I, Fig. 5g) or matched exactly the x MSPG unique value (i.e., DeDr_lT / layer I, Fig. 5m). In contrast, the statistical methods succeeded best in providing x estimates similar to those of the model of Couvreur et al. (2012) in case of a shallow groundwater table and at low T only 5 ( Fig. 5a-c and g-i), thus when water availability was high and root compensation was low. In these cases, x I_AM was included in the estimated x I_MSPG range and the mean absolute difference (MD) between x J_AM and x MSPa MFV was equal to 8.6 %.

Results and discussion
This difference was the greatest (129.2 %) for experiment DeDr_lT when HR was simulated by the analytical model ( Fig.   5m-o).
Considering eight soil layers instead of three added uncertainty in the assessment of their relative contribution to T as 10 determined by the MSPG method: the estimated probability ranges increased in most of the cases (results not shown).
However it considerably improved the results of the MSPa method: the mean absolute difference between x J_AM and the most frequent x MSPa value was equal to 4.7 % for the scenarios with a shallow groundwater table and low transpiration rate and equal to 52.1 % in case of HR (Table 3).
Independent of the number of defined soil layers, lowering the value of increment to 5 % in the MSPG method refined the 15 analysis where the probability distribution was already narrow (i.e., in the case of a well identified x MSPG value, e.g., Fig. 5m) while it produced distributions that were flatter and contained less gaps when no clear solutions had emerged before (results not shown). Artificially increasing the value of tolerance had the consequence that more solutions to Eq. (10) were found for each experiment / transpiration value / layer combinations and vice versa (results not shown). An increase or decrease of a factor 2 of the number of runs as well as the number of runs to be discarded from the analysis had only a marginal impact on 20 the density distribution curves obtained with the MSPa method in the case of three or four soil layers.
The modelling exercise illustrated the disparities of outcome between the graphical method on the one hand and the statistical and mechanistic methods on the other: there simply cannot be a single or multiple "root water uptake depths" but rather a continuous RWU profile (AM) or statistical solutions of contribution to transpiration (MSPG and MSPa). Significant changes of δ Ti do not necessarily mean important changes in the depth of RWU but rather slight (nevertheless significant) 25 modification of the RWU profile. The authors believe that the relatively novel statistical tools MSPG and MSPa presented in this review should be therefore preferred over the GI method, especially since the two former are available as user-friendly programs and packages and do not require significant computing time, therefore can be run locally on a personal computer.
As highlighted in this series of virtual experiments, the Bayesian method showed much more convincing results than the method of Phillips and Gregg (2003), especially in the case of eight soil layers, illustrating the interest of reaching the best 30 vertical resolution and maximizing the number of identified potential sources.
One can also show from this inter-comparison of methods that labelling of soil water in either 18 O or 2 H has potentials for improving the different methods presented here theoretically if water is taken up by the roots from the labeled region Biogeosciences Discuss., doi: 10.5194/bg-2016-410, 2016 Manuscript under review for journal Biogeosciences Published: 5 October 2016 c Author(s) 2016. CC-BY 3.0 License.
predominantly. However this was never the case looking at the results of the analytical model. A dual isotope ( 18 O or 2 H) labelling pulse experiment that would artificially disconnect the strong link between δ 18 O and δ 2 H would on the other hand much more constrain the inverse problem and provide accurate estimate of contribution of S to transpiration flux.

Isotopic assessment of hydraulic redistribution 5
HR (e.g., HL) has been observed using isotopic measurements in a number of studies (e.g., Caldwell and Richards, 1989;Dawson, 1993;Kurz-Besson et al., 2006). However, in contrast to nondestructive "traditional" methods allowing for direct monitoring of redistribution dynamics (i.e., psychrometry, time domain reflectometry, and frequency domain capacitance, Brooks et al., 2002;Richards and Caldwell, 1987;Wan et al., 2000), isotopic methods provide a destructive and indirect assessment. These methods are based on (i) labelling of soil or roots of deep-rooted plants at a given depth in the 10 soil or at a certain location in the experimental field and (ii) measuring the δ Ti of plants not having access to labeled water (i.e., of which the roots do not reach the isotopic labeled depth or location). When HR occurs, the xylem sap water (of measured isotopic composition δ Ti,m ) of these plants can be conceptualized as a mixture of antecedent soil water (at natural isotopic abundance) and isotopically enriched water released to the soil by the deep-rooted vegetation. From simple mass balance at the release location, δ S at a given depth z in the soil and at time (t+Δt) deviates from that at time t as a function of 15 the (negative) S (i.e., HR or HL) at time t and change of soil volumetric water content (θ): If θ and δ S at times t and t+Δt, and δ Ti are measured, the water volume transported by the roots (V HR , L 3 ) can be calculated knowing the volume of soil representative of the hydric and isotopic measurements (V, L 3 ). Note that HR is observable at a certain soil depth if and only if uptake and release locations in the soil have distinct water isotopic compositions. Finally, the 20 obtained volume can be compared with the water volume transpired by the vegetation on the following day.
To the authors' knowledge, no precise observation (other than the study of Zegada-Lizarazu and Iijima, 2004) of change of soil water isotopic composition has been attributed with certainty to hydraulic redistribution and simultaneously provided amount of water involved in the process. Such observations however should be feasible under controlled experimental conditions where (i) the initial soil water isotopic profile before labelling is known and (ii) natural isotopic changes (due to, 25 e.g., soil redistribution and moisture input from a precipitation event) can be avoided, and (iii) the lateral heterogeneity of soil water and isotopic composition profiles can be minimized (see for instance the setups of Armas et al., 2012;Querejeta et al., 2012). As highlighted in section 2.3, HR can be conceptualized as a negative S (Eq. (4b)) and should therefore be exempt of isotopic fractionation. However, to the authors' knowledge this point has not yet been proven experimentally.

High frequency isotopic data and sampling strategies
For determination of δ S , soil profiles are usually destructively sampled, typically with an auger down to a depth of a few centimeters (Rothfuss et al., 2010) to a few meters (Moreira et al., 2000) (see Table 1), depending on the depths of the root system and of the water table. The sampling depth interval should, when possible, match the exponential decrease of isotopic composition (Wang et al., 2010) due to fractionating evaporation and it should capture sudden variations with time at the soil 5 surface due to precipitation, i.e., be minimal at the surface and maximal deeper in the soil profile where isotopic dynamics are less pronounced. A minimal sampling interval at the surface is also crucial as it provides the isotopic composition of the layer contributing the most to transpiration in the case of a low T flux (e.g., morning transpiration) under non-limiting water availability. Not measuring this maximum soil isotopic composition (between precipitation event) can lead to a situation where source partitioning is not feasible from isotopic measurements. Under field conditions (i.e., ~95 % of the studies 10 reviewed in this work, summarized in Table 1) soil material is generally not a limiting factor, thus can be sampled twice or thrice to average out or characterize lateral heterogeneity without significant disturbance of the soil (Leroux et al., 1995).
Water from plant and soil materials is predominantly extracted by cryogenic vacuum distillation (Araguás-Araguás et al., 1995;Ingraham and Shadel, 1992;Koeniger et al., 2011;Orlowski et al., 2013;West et al., 2006). This consists in (i) introducing the plant or soil sample into an extraction flask attached to one end of the extraction line, while at the other end a 15 collection tube is connected, (ii) freezing the sample by immersing the collection flask into liquid nitrogen (temperature ~ -200°C), (iii) pumping the extraction line down to a pressure of ~10 -3 mbar, (iv) heating the sample to a certain temperature ( ~60 < T < ~100°C) depending on its nature while immersing the trap into liquid nitrogen. The water vapor produced condenses in the trap following a stepwise procedure (~lasting one to a few hours), in order to avoid condensation elsewhere on the water vapor path between sample and collection trap. Accuracy of this extraction method was shown to be maximal at 20 higher water content and for sandy soils and lower for soils with high clay content. In the latter case, extraction times should be longer and temperatures higher to mobilize water strongly bound to clay particles, which has a distinct isotopic composition from that of pore "bulk" water (Araguás-Araguás et al., 1995;Ingraham and Shadel, 1992;Oerter et al., 2014;Sofer and Gat, 1972). In other studies, plant and soil waters are extracted following azeotropic distillation with kerosene as solvent (e.g., Brunel et al., 1995;Thorburn and Ehleringer, 1995), or direct equilibration with CO 2 (Asbjornsen et al., 2007) 25 following the method of Scrimgeour (1995), or else the mild vacuum method (Dawson and Pate, 1996;Jeschke and Pate, 1995).
Certainly one of the main limitations of all isotopic approaches for quantifying RWU and HR is the destructive character of isotopic sampling (see section 3.1) and associated offline analyses (sections 2.2 and 2.3). This usually leads to poor spatial (maximum a few cm 2 ) as well as temporal (minimum hourly) resolution of the inferred results, when comparing with 30 measuring frequency of other soil and plant state variables, e.g., soil water content and potential, and leaf water potential very different water residence time is mixed. Similarly, given the expected high lateral and temporal variability of the HR process, the representativeness of δ S should be questioned for soils, in particular when combined with 1D models.
Recently developed methods take advantage of laser-based spectroscopy which allows on-line and continuous isotopic measurements in the gas phase. These methods rely on coupling a laser spectrometer with specific soil gas sampling probes consisting of gas-permeable microporous polypropylene membranes or tubing. These membranes or tubing exhibit strong 5 hydrophobic properties, while their microporous structures allow the intrusion and collection of soil water vapor. Several authors (Gaj et al., 2015;Gangi et al., 2015;Herbstritt et al., 2012;Rothfuss et al., 2013;Sprenger et al., 2015;Volkmann and Weiler, 2014) could determine the soil liquid water isotopic composition in a nondestructive (yet invasive) manner from that measured in the collected soil water vapor considering thermodynamic equilibrium between vapor and liquid phase in the soil. In contrast to "traditional" isotopic methods, these novel isotopic monitoring methods have also the distinct 10 advantage of determining soil liquid water isotopic composition at very low water content, since water vapor, in contrast to soil liquid water, is not limiting for analysis. These novel methods allow a vertical resolution down to 1 cm and an approximately hourly time-resolution. However, they do not allow horizontal resolution along the tube and are greatly sensitive to the carrier gas used (Gralher et al., 2016). In their opinion papers, McDonnell (2014) and Orlowski et al. (2016) also urged for a comparison between methods, which was addressed by Gaj et al. (2015) and Pratt et al. (2015). 15 Leaf and plant gas chamber systems provide indirect means for a nondestructive determination of δ Ti , i.e., by either assuming full steady-state conditions at the evaporative sites of the leaves (δ Ti = δ T ) (e.g., Dubbert et al., 2014;Volkmann et al., 2016b). In the coming years, effort should be made towards developing novel methods for a direct and nondestructive determination of δ Ti based on the use of gas-permeable membranes, which was recently initiated for trees (Volkmann et al., 2016a). This should be further investigated to test applicability to other (non-woody) plant species. This will imply the major 20 challenge of not disrupting the water columns in the active xylem vessels when installing such a membrane-based system.
Another potential issue to be investigated is the species-specific extent of water exchange between xylem and phloem conductive tissues which might lead to isotopic "contamination" of the xylem sap water.

Call for a coupled experiments-modelling approach for determination of plant water sources and redistribution on the basis of isotopic data 25
In order to fully benefit from the potential of water stable isotopologue analysis as tools for partitioning transpiration flux, the authors call for a generalization of coupled approaches based on the confrontation of experimental data with a physically based understanding of RWU processes.
Simple analytical models, such as the formulation of Couvreur et al. (2012), can be applied and confronted with isotopic data. In comparison with statistical tools, such physical models provide profiles with high spatial resolution and lower 30 uncertainty, on the condition that all required (isotopic) data is available. We recognize that in comparison with the statistical and conceptual methodologies presented in this review, using a physical (analytical or numerical) model implies the Biogeosciences Discuss., doi: 10.5194/bg-2016-410, 2016 Manuscript under review for journal Biogeosciences Published: 5 October 2016 c Author(s) 2016. CC-BY 3.0 License. measurements of additional state variables to be fed as input to the model, and of one parameter (K plant ) (when considering the assumption K plant = K comp valid, see Appendix B). Some of these variables are laborious to obtain (e.g., RLD) or not straightforward to measure (H S , H L , and T)especially in the fieldbut are mandatory to be able to determine contributions to T across a set of identified water sources. In addition, they are necessary to gain insights into soil-plant interactions, e.g., dynamics of root function (active versus non-active roots in the soil profile) in water uptake and thus quantify the 5 disconnection between measured RLD and the prognostic variable SSF (see Appendix B1) For doing this, controlled conditions in state-of-the-art climatic chambers are ideal, as they allow reducing the inherent spatial heterogeneity present under natural conditions and, thus, the deconvolution of environmental effects on RWU. Experimental facilities that not only control atmospheric forcing (soil upper boundary conditions for latent and heat flow), but impose lower boundaries for the soil compartment (e.g., drainage and capillary rise dynamics) and provide means to close the hydrological balance are 10 required. Moreover, macrocosm experiments (~m 3 scale) should be favored over mesocosm (~dm 3 scale) experiments to avoid or reduce inherent side effects that would ultimately hamper mimicking natural conditions. Biogeosciences Discuss., doi: 10.5194/bg-2016-410, 2016 Manuscript under review for journal Biogeosciences Published: 5 October 2016 c Author(s) 2016. CC-BY 3.0 License.

CONCLUSION
Root water uptake is a key process in the global water cycle. More than 50% of total terrestrial evapotranspiration crosses plant roots to go back to the atmosphere (Jasechko et al., 2013). Despite its importance, quantification of root water uptake remains difficult due to the opaque nature of the soil and the spatial and temporal variability of the uptake process.
Water stable isotopic analysis is powerful and valuable tool for the assessment of plant water sources and for the 5 identification of hydraulic redistribution. In an inverse modelling framework, isotopic analysis of plant tissues and soil also allow for obtaining species-specific parametrization of physically-based analytical and numerical RWU models. They provide at the plant scale a unique way to tackle the difficulty of disentangling actual RWU profiles with root traits and characteristics.
In this review we tried to highlight the importance of systematically reporting uncertainties along with estimates of 10 contribution to T of given plant water sources. The inter-comparison exercise could quantify the impact of the definition of the plant water sources (i.e., whether they are spatially disjoint or not and whether their isotopic compositions values are significantly different or not) on the outcome of the two end-member mixing model. The inter-comparison also illustrated the limitations of the graphical inference method and the multi-source mixing model of Phillips and Gregg (2003), whereas it underlined the performance of the Bayesian approach of (Parnell et al., 2013), which uses a more rigorous statistical 15 framework, if the number of considered water sources matches the number of isotopic measurements in the soil profile.

8
The color-shaded areas depict the results of 1000 model runs where for each input data variable (soil water potential, 9 δ 18 O, and root length density -RLD) a single offset randomly selected between -5 and +5 cm, -0.2 and +0.2 ‰, and -10 0.1 and +0.1 cm cm -3 respectively for each variable was added to the initial values reported in Table 2   where the soil has a shallow (deep) groundwater table while "Dr" and "We" stand for when the soil is dry or wet at the surface (e.g., shortly after a rain event). Suffices "lT" and "hT" refer to "low" and "high" transpiration rate simulations. "*" refers to when hydraulic redistribution is simulated by the analytical model, leading to a negative x. Error bars 10 refer to either one standard deviation (for the RWU analytical model) or one standard error (for the TM approach).   (Table 1) and the δ Ti simulated by the model of Couvreur et al. (2012) (i.e., δ Ti_AM, given above each plot along with its standard deviation). Tolerance of the MSPG was set equal to the standard deviation of δ Ti_AM . "Sh" ("De") stands for the virtual experiments where the soil has a shallow (deep) groundwater table while "Dr" and "We" stand for when the soil is dry or wet at the surface (e.g., shortly after a rain event). Suffices "lT" and "hT" refer to "low" and "high" transpiration rate simulations. The colored vertical lines give x I_AM , x II_AM , and  Vienna-Standard Mean Ocean Water (V-SMOW) hydrogen or oxygen stable isotopic ratio Root water uptake sink term, Root water uptake sink term under uniform soil water potential distribution, compensatory root water uptake sink term, total root water uptake sink term as simulated by the analytical model of Couvreur et al. (2012) Standard sink fraction Time, time step Transpiration flux Contributive proportion to transpiration, source j contributive proportion to transpiration, continuous and integrated (layer J) contributive proportions to transpiration as simulated by the analytical model of Couvreur et al. (2012), integrated (layer J) contributive proportions to transpiration as determined by the statistical approaches of Phillips and Gregg (2003) and Parnell et al. (2010). Contributive proportion to transpiration under conditions of uniform soil water potential Soil depth, soil depth of layers j and j+1, thickness of soil layer j, depth of the root system, "mean root water uptake depth"   O, δl, δv, δl-v, δsource, δsurf, δsim, δS, δS,j, δS,J, δA, δB, δTi, δTi,m, δTi_AM, δE, δT Equilibrium and kinetic isotopic fractionation factors, hydrogen and oxygen equilibrium isotopic fractionation factors, hydrogen and oxygen kinetic isotopic fractionation factors Water stable isotopic composition, water hydrogen and oxygen stable isotopic compositions, liquid, vapor, liquid-vapor interface, source, soil surface, and simulated water isotopic compositions, soil water isotopic composition, soil layer j and J water isotopic composition, sources A and B water stable isotopic compositions, isotopic composition of xylem sap water at the plant tiller, isotopic composition of xylem sap water measured at the plant tiller, isotopic composition of xylem sap water at the plant tiller as simulated by the model of Couvreur et al. (2012), isotopic composition of transpiration Residual error term Soil volumetric water content Volumetric mass of water Sandard errors associated with the measurements of x, δA, δB, δTi and estimated uncertainty of δTi_AM as simulated by the analytical model of Couvreur et al. (2012), error associated with the estimation of the contributive proportion to T of water source A in the case of two distinct sources S can become positive, and water is released to the soil (i.e., hydraulic redistribution -HR occurs). From Eq. (B4), it can be concluded that HR will preferably occur when q Ti is small and when large soil water potential gradients exist. Plant root hydraulic characteristics will control compensation through the K comp term. The importance of the compensatory RWU term has been discussed in the literature for a long time (e.g., Jarvis, 1989). Except if plants activate specific mechanisms to avoid it, compensation always takes place under natural conditions due to the spatially heterogeneous distribution of soil water 5 potential (Javaux et al., 2013).
A simplifying hypothesis that can be made (Couvreur et al., 2014;Couvreur et al., 2012) is to consider that K plant and K comp are equal, which substituted in Eq. (B4) leads to: Finally, the uptake of water stable isotopologues, i.e., the "isotopic sink term" (S i [M T -1 ]) is defined as: 10 where C [M L -3 ] is the water isotopic concentration.

B2: Running the model for the inter-comparison
The root water uptake (S AM ) depth profiles and corresponding δ Ti_AM were simulated using the model of Couvreur et al. from each of the integrated soil layers J (J ≤ III or J ≤ VIII) were calculated. In order to account for uncertainty of the input data (i.e., total soil water potential and oxygen isotopic 20 composition H S and δ S , and root length density RLD), the model was run a 1000 times where a single offset randomly selected between -5 and +5 cm, -0.2 and +0.2 ‰, and -0.1 and +0.1 cm cm -3 was added to the initial values (reported Table   2) of H S , δ S , and RLD, respectively. By doing this we obtained a posteriori distributions of S AM and corresponding δ Ti_AM