Reviews and syntheses: Isotopic approaches to quantify root water uptake: a review and comparison of methods

. Plant root water uptake (RWU) has been docu-mented for the past ﬁve decades from water stable isotopic analysis. By comparing the (hydrogen or oxygen) stable isotopic compositions 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), studies were able to determine the relative contributions of these water sources to RWU. In this paper, the different methods used for locat-ing/quantifying relative contributions of water sources to RWU (i.e., graphical inference, statistical (e.g., Bayesian) multi-source linear mixing models) are reviewed with em-phasis 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 while it underlines

Abstract. Plant root water uptake (RWU) has been documented for the past five decades from water stable isotopic analysis. By comparing the (hydrogen or oxygen) stable isotopic compositions 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), studies were able to determine the relative contributions of these water sources to RWU.
In this paper, the different methods used for locating/quantifying relative contributions of water sources to RWU (i.e., 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 while it underlines the performance of one Bayesian mixing model. 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, the authors call for a development of approaches coupling physically based RWU models with controlled condition experimental setups.

Introduction
Understanding how the distribution of soil water and root hydraulic architecture impact root water uptake (RWU) location and magnitude is important for better managing plant irrigation, developing new plant genotypes more tolerant to drought or tackling ecological questions in water-limited ecosystems, such as the competition for soil water by different plants (Javaux et al., 2013).
RWU -defined as the amount of water abstracted by a root system from soil over a certain period of time -is principally driven by transpiration flux taking place in the leaves. 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 stress hormonal signals are related to the soil water potential distribution and to the plant hydraulic architecture (Huber et al., 2015). The distribution of RWU is very variable in time and space, and depends on the presence of roots and their 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, i.e., the ability of the soil to provide water at the plant imposed rate (Couvreur et al., 2014): 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 rhizosphere, at the soil-root interface and in the root system (Steudle and Peterson, 1998). The actual RWU Published by Copernicus Publications on behalf of the European Geosciences Union. R L D �cm cm -3 � Figure 1. Some examples of root water uptake sink term (S, in d −1 ) profiles (blue lines) conceptualized as the sum of S uniH (green lines), the root water uptake term proportional to root length density RLD (dotted black line) and the compensatory root water uptake (S comp , red lines). (a) S comp = 0 (no root compensation) leading to S = S uniH . (b) S comp = 0 and becomes negative towards the surface but remains smaller (in absolute terms) than S uniH . (c) S comp = 0 and is negative at the surface, while it is greater than S uniH for z > −0.08 m. In the last case, S is negative at the surface, meaning hydraulic lift is observed.
profile is thus a combination of different aspects: the root's ability to extract water (characterized by the number 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, e.g., adaptive root growth, adaptive root conductivity or 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 length densities (RLD (L L −3 ), usually expressed in cm root per cm 3 soil), 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 T −1 L −3 ), described as a sink term in the Richards equation. According to Couvreur et al. (2012), root compensation is defined 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: S(x, y, z) = S uniH (x, y, z) + S comp (x, y, z), where x, y and z are the three-dimensional (3-D) spatial coordinates, S uniH is a term proportional to the root distribution and S comp is the 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 1-D space. When there is no compensation (S comp (x, y, z) = 0), the RWU distribution follows the root distribution (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 (called "Hydraulic redistribution" or "Hydraulic lift" in this particular case, Caldwell and Richards, 1989;Dawson, 1993;Kurz-Besson et al., 2006). Despite its importance, datasets with measurements of RWU are lacking. This is 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 in 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 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 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 al., 2000). More recently, Zarebanadkouki et al. (2012) were able to measure for the first time RWU in porous media by combining a tracer experiment (i.e., deuterated water) monitored by neutron tomography with inverse modeling of a transport equation. Yet this was performed under controlled conditions, while there is no standard method to monitor three-dimensional water  Nadezhdina et al. (2010Nadezhdina et al. ( , 2012Nadezhdina et al. ( , 2015 used sap flow measurements in roots to quantify hydraulic redistribution. 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 isotopes, water stable isotopologues ( 1 H 2 H 16 O and 1 H 2 18 O) have been frequently used to identify and quantify RWU in soils through the measurements of their natural (and artificial) abundances. Methods include simple graphical inference to more sophisticated statistical methods, i.e., twoend members and multi-source linear mixing models. While the former attempts to locate the "mean root water uptake depth" 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 has three objectives: (i) performing a literature review on the use of water stable isotopes to assess RWU; (ii) presenting the methods for translating the isotopic information into RWU profiles (i.e., graphical inference and statistical multi-source linear mixing models); and (iii) comparing these methods with a series of virtual experiments differing in the water and isotopic statuses in the soil and the plant. Prior to the review and inter-comparison, 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 evoke opportunities offered by novel isotopic monitoring tools which provide unprecedented high-frequency isotopic measurements, and call for a development of approaches making use of physically based models for RWU determination.

Flow of isotopologues in the soil-plant system
In a study that laid the basis for future work in isotopic ecohydrology, Zimmermann et al. (1967) provided a steadystate analytical solution for soil water isotopic composition (δ S , expressed in ‰ relative to the Vienna Standard Median Ocean Water international (VSMOW) isotope reference scale, Gonfiantini, 1978) 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 (E) 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. Furthermore, by conservation of mass, the isotopic composition of evaporation equals that of its source (e.g., groundwater), i.e., δ E = δ source . A profile is obtained (Fig. 2a, dark blue line) whose exponential shape depends on boundary conditions, i.e., the source water and 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. Barnes and Allison (1983) extended this formulation to a non-saturated sand column evaporating at isotopic steady state. 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. In the "vapor region", relative humidity generally is still close to unity for sand total water potential below 15 bars. 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" 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, light blue line). Finally, note that Rothfuss et al. (2015) argued that, at transient state (δ E = δ source ), the maximal isotopic enrichment in the soil profile might not point to the location of the evaporation front. Instead, they proposed that the depth where the steepest gradient in the isotopic profile is observed 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 greater than two and lower than six, 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 evaporation front (2 < δ 18 O-δ 2 H slope < 6) (Sprenger et al., 2016). 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. Finally, under natural conditions, the δ S profile is not solely a result of isotopic fractionation, but is also highly impacted both spatially and temporally by input precipitation isotopic composition through modification of the upper boundary condition (δ surf ).
As opposed to the removal of water vapor by evaporation, RWU has been described in a number of studies and over a wide 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  Figure 2. (a) Simulated soil water isotopic composition (δ S ) profiles for a water saturated (dark blue line) and unsaturated (light blue dotted line) soil following Barnes and Allison (1983). Indices "surf" and "EF" refer to soil surface and evaporation front. "vapor" and "liquid" regions refer to soil regions where water flow occurs predominantly in the liquid and vapor phases, respectively. (b) Illustration of the graphical inference method for determining the mean root water uptake depth (z). "Ti" stands for the sap xylem water at the plant tiller. Case 1: one unique solution is found; case 2: more than one solution is found. A smoothed profile is designated by the symbols. Thez range is indicated by the gray horizontal stripes. composition due to evaporation exists, the root system takes up water at different depths, thus having different isotopic compositions.
Assuming that water transport time in roots is negligible, the isotopic concentration of the xylem sap water at the root tiller (C Ti (M L −3 )) can be modeled as the weighted average of the product of the soil water isotopic concentration (C S (M L −3 )) and S(x, y, z): with J Ti (L 3 T −1 ) the xylem sap flux at the root tiller. Following Braud et al. (2005), with δ S (x, y, z) (-, expressed in ‰) the isotopic compositions of soil water at coordinates (x, y, z). Mostly, a onedimensional description of RWU 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): 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 (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 Biogeosciences, 14, 2199Biogeosciences, 14, -2224Biogeosciences, 14, , 2017 www.biogeosciences.net/14/2199/2017/ 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 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 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 O 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. Finally, plant water samples will, similarly to soil water samples, fall onto an "evaporation line" of a slope lower than eight in a two-dimensional (δ 18 O, δ 2 H) space (Javaux et al., 2016).

Literature review
By entering the search terms (("root water uptake" or ("water source" and root) or "water uptake") and isotop * ) into the ISI Web of Science search engine (www.webofknowledge.com), 159 studies published in the last 32 years were identified (see a listing of all studies in the Supplement). Cumulative number of articles as a function of publication year follows an exponential shape: on average over the period 1985-2014, number of publications per year increased for about 0.3 and reached 8 (2014). In both years 2015 and 2016, the isotopic method for locating or partitioning water sources to RWU gained significantly more attention with 20 publications per year (Fig. 3a).
When sorting plant species simply by their form and height, it appears that trees are the most studied group of plants (present in about 60 % of the studies), followed by annual and perennial grasses (21 %) and shrubs (e.g., desert and mangrove species, 21 %) (Fig. 3b). Only 15 % of the publications study RWU in agricultural systems (e.g., maize, wheat, millet, rice), which is reflected by the small portion of peer-reviewed journals of which the category is listed under "Agronomy and Crop Science" (8 %) by Scimago Journal & Country Rank (www.scimagojr.com). This is a rather surprising finding given the fact that drought stress is considered as a major threat for crop yields and that RWU is a crucial mechanism to sustain drought periods. "Soil Science" is a relatively underrepresented category with 8 % as well. This is corroborated by the fact that 27 % of the studies do not report any information about soils (e.g., texture, FAO class, structure, particle size distribution, or physical properties) (Fig. 3c). In comparison, the "Ecology" and "Ecology, Evolution, Behavior and Systematics" categories are significantly more represented with 22 % of the studies altogether.
Four classes of methods for RWU analysis on basis of isotopic information emerged from our analysis (Fig. 3d). In a first one, representing 46 % of the studies, RWU is either located in a specific soil layer using the method of "direct inference" (Brunel et al., 1995) or in some water pool (or water "source", not to be mistaken with the concept of water source defined in the previous section), e.g., groundwater, soil water, or rainwater (Andrade et al., 2005;Beyer et al., 2016;Roupsard et al., 1999). In a second class (32 % of the studies), relative contributions of at least three water sources to RWU are determined using multi-source mixing models (e.g., IsoSource, Phillips and Gregg, 2003, representing 21 % of the studies; SIAR, Parnell et al., 2013, 5 %;MixSir and MixSIAR, Moore and Semmens, 2008, 2 %). In a third class (18 % of the studies), relative contributions of two particular water sources (e.g., water in two distinct soil layers, or groundwater versus recent precipitation) to RWU are calculated "by hand" with a two-end-member linear mixing model (Araki and Iijima, 2005;Dawson and Pate, 1996;Schwendenmann et al., 2015). Note that classes two and three (representing 50 % of the studies) are both based on end-member mixing analysis (EMMA) (Barthold et al., 2011;Christophersen and Hooper, 1992) and will be further pooled into "statistical approach" in Sect. 3.2 of this study.
Note that all methods have in common the use of an inverse modeling approach: the RWU distribution is obtained by optimizing model input parameters until the simulated δ Ti and/or the simulated δ S profiles fit to the isotopic measurements. One important feature of the three first classes of methods is that they consider soil water isotopic transport flow to be negligible for the duration of the experiment. Numerical models such as HYDRUS-1D and SiSPAT-Isotope on the other hand take this into account in the computation of RWU profiles. The first three methods also differ from the last one by the fact that they only give fractions of RWU instead of absolute RWU rates changing in time and space.
Sixty-one percent of the studies based their estimation of location or quantification of relative contributions on measurement of either δ 2 H or δ 18 O, i.e., in a single isotope framework, while 25 % used both δ 2 H and δ 18 O (i.e., in a dual isotope framework). In the remaining studies, both isotopic compositions were measured and used to provide two separate estimates of relative contribution distributions even though δ 2 H or δ 18 O distributions were strongly linked (see Sect. 2). This last approach is in the present study referred as "double single" (see the Supplement). The vast majority of the studies (82 %) took advantage of natural isotopic abundances, while the rest (18 %) applied labeling pulses to the soil (either in the profile or at the soil boundaries, e.g., the soil surface and groundwater) to infer RWU from uptake of labeled water.
To summarize, we observe that isotopic analyses have mainly been used up to now to assess water sources under natural ecosystems mainly using statistical approaches. On the opposite, these techniques have not been used much to investigate RWU of crops. It is also observed that the use of water isotope composition datasets combined with explicit physical models is lacking. In the next sections, we analyze the main methods currently use to retrieve RWU with water isotopic compositions and compare the different methods. Table 1 summarizes 21 particular isotopic studies that use either one of the first three classes of methods (i.e., accounting for about 96 % of the published studies), while class four (physically based RWU models) will be treated separately in Sect. 4 of this study. These 21 studies were chosen according to either the number of citations and contribution importance (for studies published before 2015) or to the novelty of the publications (publication year ≥ 2015).

Graphical inference (GI)
This straightforward approach first proposed by Brunel et al. (1995) and applied by, e.g., Leroux et al. (1995), Weltzin and McPherson (1997) (Table 1) and elsewhere by Midwood et al. (1998), Armas et al. (2012), and Isaac et al. (2014) (see the Supplement) defines the "mean root water uptake depth" z as the depth where δ S = δ Ti .z conceptually indicates the soil depth where the plant root system, represented as one unique root, would extract water from.
There are cases wherez cannot be unambiguously identified (e.g.,z 1 andz 2 of case 2, Fig. 2b) due to the nonmonotonic 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 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 in Fig. 2b where a uniform θ profile is assumed, the δ S,J profile intersects with the vertical line of value δ Ti deeper than for the  www.biogeosciences.net/14/2199/2017/ initially non-monotonic δ S profile, i.e.,z (case 2, integrated δ S profile) <z 2 <z 1 . Some authors rule out solutions in the case of multiple mean RWU depths, e.g., by excluding thē z solutions where soil water content is low and/or soil water potential is 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 equivalent to the isotopic composition "sensed by the plant" only if the root profile is homogeneous, i.e., when RLD is constant over depth in that particular soil layer J .
The graphical inference method may not only providez 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 the uncertainty in determiningz is. Figure 2b illustrates this with estimated minimum and maximumz for the monotonic δ S profile and for the vertically averaged profile. In the latter case, the possible range ofz is the largest. These ranges give first quantitative indication of variance aroundz. Finally, for a complete "graphical assessment" of the variance ofz, one should consider the uncertainty associated with measurements of the δ S profile as well (not shown here; for a complete assessment of errors associated with determination of δ S , see Sprenger et al., 2015).

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 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) 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 with J A , J B , and J Ti (L 3 T −1 ) (respectively, J i A , J i B , and J i Ti (M T −1 )), and the fluxes of water (or isotopologues) originating from water sources A and B, and at the plant tiller. C A , C B , and C Ti (M L −3 ) are water sources A and B, and xylem sap water measured isotopic concentrations. By introducing x = J A /J Ti and following Eq. (3), Eq. (6b) becomes 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., x and (1−x), respectively. The error associated with the estimation of x (σ x (-, expressed in ‰)) can be calculated following Phillips and Gregg (2001): with σ δ A , σ δ B , and σ δ Ti the standard errors associated with the measurements of δ A , δ B , and δ Ti , respectively. The sensitivity of Eq. (8b) to different values of σ δ A , σ δ B , and σ δ Ti can be tested by considering either minimal possible errors, i.e., the analytical precision of the isotopic analyzer (e.g., isotope ratio mass spectrometer, laser-based spectrometer), or by taking additional errors involved with sampling procedure and vacuum distillation technique (see, e.g., Rothfuss et al., 2010) into account. Equation (8b) also shows that, independently of the values considered for σ δ A , σ δ B , or σ δ Ti , σ x is inversely proportional to 1/(δ A − δ B ), indicating that the two end-members should have the greatest possible isotopic dissimilarities 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 compositions, and that the difference is considerably larger than the precision of the isotopic measurements. Figure 4 shows for example that when (i) x is evaluated at 10 % and (ii) σ δ A , σ δ B , and σ δ Ti are estimated as equal to 0.02 ‰ (dark blue solid line), (δ A − δ B ) should be greater than 0.75 ‰ (in absolute terms) in order to reach a σ x value lower than 5 %, i.e., more than 37 times the error made in δ A , δ B , and δ Ti . To obtain the same standard error for x in the case of a higher standard error in the estimation of δ A , δ B , and δ Ti (e.g., σ δ A , σ δ B , and σ δ Ti = 0.1 ‰), (δ A − δ B ) should be greater than 3.00 ‰ (in absolute terms). This difference becomes much greater for σ δ A , σ δ B , and σ δ Ti = 1.00 ‰ and reaches 42 ‰ (not shown in Fig. 4). This certainly highlights the advantage of artificially labeling 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 portions 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 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 were able to find evidence against the 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 4 also illustrates that for given (δ A − δ B ), σ δ A , σ δ B , and σ δ Ti values, the "optimal x value" for a low σ x is 50 % (shown by the orange lines). Table 1 displays a sample of studies that used the twoend-member mixing approach. Authors were able to 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) were for instance able to locate the dominant source for RWU, i.e., groundwater for their mountain and a floodplain test site and water from the soil between 0.3 and 0.4 m depths 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, or else calculated the ratio of geometrical distances between δ Ti and δ A and between δ Ti and δ B in dual isotope (δ 18 O and δ 2 H) space (Bijoor et al., 2012;Feikema et al., 2010;Gaines et al., 2016). 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 to as "D" in Table 1) will or should gain in importance in isotopic studies. This is especially useful when (i) under natural conditions the δ 18 Oδ 2 H slope is not constant over depth (Sprenger et al., 2016) or (ii) in the context of pulse labeling experiments, which can artificially change the value of the δ 18 O-δ 2 H slope at given locations in the soil profile. In these cases, two independent mixing equations are obtained, 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 the soil profile, Eq. (7) becomes with N the number of plant water sources (e.g., soil layers) and N j =1 x j = 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, based on background information or knowledge, some of these solutions are not likely or possible. 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 number of combinations depends on the value of the contributing 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 the 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, which therefore should be carefully chosen by the user; e.g., a smaller i refines the analysis. For this, the IsoSource program (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 nonuniqueness problem and provided diffuse patterns of frequency that were difficult to interpret in some cases (e.g., in the case of a non-monotonic isotopic profile), it had the advantage over the former method of providing a systematic and quantitative assessment of ranges of relative contributions. Romero-Saltos et al. (2005) extended the model of Phillips and Gregg (2003) by constraining RWU to follow a normal distribution within a delimited 50 cm soil vertical segment of 1 cm vertical resolution and centered aroundz, the mean RWU depth. The location of this section and thusz is also obtained by mass balance from inverse modeling similarly to IsoSource (see applications of Grossiord et al., 2014;Rossatto et al., 2013;Stahl et al., 2013). 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 to Biogeosciences, 14, 2199-2224, 2017 www.biogeosciences.net/14/2199/2017/ (ii) provide a optimal solution rather than ranges of feasible solutions. For doing this, they use a Bayesian framework (for details see Erhardt 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 ): x j · δ S,j + ε j .
Note that the terms of (i) trophic enrichment factor (TEF (-, expressed in ‰); see, e.g., the meta-analysis of Vanderklift and Ponsard, 2003) and (ii) isotope concentration dependency Phillips and Koch, 2002) originally incorporated into 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.r-project.org/src/ contrib/siar_4.2.tar.gz) in which the initial (a priori) x j distribution is by default the Dirichlet distribution, of 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-Hastings (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 inz and change in the 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 inz was non-existent or not 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 profile time series with time series of δ Ti (indirectly obtained from the isotopic measurement of the transpired water and assuming isotopic steady state, i.e., δ Ti = δ T ) following a labeling pulse (see Table 1 for details on the study).

Inter-comparison of methods
We tested and compared the different methods (GI, TM, MS) during a series of virtual experiments in a single isotope framework (δ 18 O) and at natural isotopic abundance. Mean RWU depths (provided by the GI method) and x j distribution (provided by the two-end-member and multi-source mixing models) were determined from the δ S profile and the δ Ti value. For each virtual experiment the δ S profile was prescribed to the different methods, while δ Ti was calculated with the physically based analytical RWU model (referred to as "Couv") of Couvreur et al. (2012).
It has been proved (Couvreur et al., 2012) that this model gives similar results to a 3-D physically based model with detailed descriptions of the root architecture and of the water flow in soil and roots. In that sense, this is the best current model existing nowadays to simulate water fluxes in a soilplant system (based on biophysical considerations). Other current models make assumptions or use empirical relations to predict RWU, which are not based on bio-physical considerations only (Jarvis, 2011;Simunek and Hopmans, 2009). Obviously, we do not mean that the model of Couvreur et al. (2012) gives the reality, but rather the best estimate of the water flow based on our physical knowledge.
The inter-comparison of models was performed using a single isotope ( 18 O) approach as the focus here was the differences of outcomes rather than the impact of the input isotopic data on these results. The reader is referred to Appendix B1 for a description of the model of Couvreur et al. (2012) and to Appendix B2 on how it was implemented for the inter-comparison.

Scenario definition
We developed eight virtual plausible scenarios of soil-plant systems under different environmental conditions. For each scenario, we set one total soil water potential (H S ) profile and one soil water oxygen isotopic composition (δ S ) profile. These profiles resulted from the combination of a lower boundary condition, i.e., the depth of the groundwater table, and an upper boundary condition, i.e., the soil surface water status. The groundwater table (of water isotopic composition equal to −7 ‰) 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"). For instance, for scenario "ShDr", we set the δ S profile to be maximal at the surface, due to evaporation, and minimal from −0.5 m downward, due to the shallow groundwater table location. For scenario "DeWe", on the other hand, the increase in δ S towards the surface was not monotonic due to a recent precipitation event (of water isotopic composition equal to −7 ‰). Finally, we tested two different values of plant transpiration rate (T ) and leaf water potential (H L ) with each of these four combinations (i.e., ShDr, ShWe, DeDr, and DeWe). 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"). All eight scenarios relied on a common measured root length density vertical distribution of Festuca arundinacea. Table 2 reports the input data. Note that, as hypothesized in Eq. (4b), www.biogeosciences.net/14/2199/2017/ Biogeosciences, 14, 2199-2224, 2017 Table 2. Soil, plant, and isotopic synthetic input data for the different modeling approaches (depth (z) profiles of soil water content θ , total soil water potential H S , soil water oxygen isotopic composition δ S , root length density RLD, transpiration rate T , and leaf water potential H L ) "collected" during eight virtual experiments differing in the depth of the groundwater   transpiration and sap flow rates (i.e., per unit of surface area (L T −1 )) were considered equal. The objective was not to use an advanced numerical model such as, e.g., SiSPAT-Isotope or Soil-litter-iso, to produce these scenarios, but rather use synthetic information based on (i) experimental data and (ii) expert-knowledge which would ideally illustrate the performances or limitations of the different methods.

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 layers spreading from 0 to 0.225 and from 0.225 to 2.00 m and (ii) two disjoint soil layers spreading from 0 to 0.225 and from 1.75 to 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.
For the multi-source mixing models of Phillips and Gregg (2003) (IsoSource) and Parnell et al. (2010) (SIAR), 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. IsoSource and SIAR were tested for eight soil layers (i.e., as many layers as measurement 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) as well. Increment and tolerance in IsoSource were fixed to 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 SIAR model, uncertainty associated with δ S measurements was set to 0.2 ‰ and the number of iterations was fixed to 500 000 and number of iterations to be discarded to 50 000.
For a detailed description of the inter-comparison methodology, the reader is referred to Appendix C.

Results and discussion
Figure 5 displays x Couv , the S(z)dz T /( x· y) ratios (solid colored lines) simulated by the analytical model of Couvreur et al. (2012) for the eight scenarios together with uncertainty (shaded areas) and the corresponding δ Ti_Couv (±1 SD) (for a description on how uncertainty was assessed, refer to Appendix C). In general, at high T the compensation was negligible and the S profile was mainly proportional to the RLD profile (Fig. 5b, d, f, and h). The only exception was a soil with a deep groundwater table and a dry surface, where this dry layer limited RWU (DeDr_hT). At lower transpiration demand, the S profile predicted by the Couvreur et al. (2012) model generally differed from the RLD profile (Fig. 5a, c, e, and g). This is due to the fact that the second term of Eq.  Couvreur et al. (2012) for experiments "ShDr" (soil with a shallow groundwater table and a relatively dry soil surface), "ShWe" (soil with a shallow groundwater table after a rainfall event), "DeDr" (soil with a deep groundwater table with a relatively dry soil surface), and "DeWe" (soil with a deep groundwater table and a wet soil surface). Suffices "lT" and "hT" refer to "low" and "high" transpiration rate simulations. The color-shaded areas depict the uncertainty of the model estimates on account of the precision of the measurements. The horizontal gray-shaded areas delimit the mean root water uptake layer(s) obtained by the graphical inference method. At the bottom right corner of each plot is a detail presented for z ≥ −0.10 m. Finally, results from the first term of the model of Couvreur et al. (2012) which considers uptake proportional to root length density are plotted as a dashed brown line for comparison. Note that negative x Couv means hydraulic redistribution by the roots.
tionally 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., hydraulic redistribution) was observed only for the soil with the deep groundwater table and dry upper layer (DeDr_hT, Fig. 5e). From the graphical method GI, either a single or two distinct solutions forz (displayed as gray-shaded horizontal stripes) were retrieved, depending on the monotonic/non-monotonic character of the δ S profile, and ranged between −0.02 and −0.95 m. Figure 6a displays the relative contribution to T of the uppermost layer 0-0.225 m in the case of two conjoint soil layers as computed with the TM approach and a comparison with the results of the analytical model. Except for the very last two virtual experiments with the deep groundwater table and wet upper layer at both low and high transpiration rates (i.e., DeWe_lT and DeWe_hT), there was a very good agreement between the analytical model and the two-end-member mixing model: the absolute difference between x Couv and x TM ranged between 1.5 % (ShDr_lT) and 6.3 % (ShDr_hT). For the experiment with the deep groundwater table and dry upper layer at a low transpiration rate (DeDr_lT), the TM approach estimated that x was equal to 12.3 %, while the analytical model simulated hydraulic redistribution, i.e., excluded the 0-0.225 m layer as a potential source. The significant difference between results of the two models during experiments DeWe_lT and DeWe_hT and the higher standard error associated with x TM (σ x , displayed in the form of error i.e., information on soil water isotopic composition is lacking between 0.225 and 1.75 m. "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. Suffixes "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 are either 1 standard deviation (for the RWU analytical model) or 1 standard error (for the TM approach). bars in Fig. 6) 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 Sect. 3.2.1. Figure 6b gives the relative contribution to T of the 1.75-2.00 m layer in the case of two disjoint soil layers, i.e., when not all potential water sources are accounted for in the calculation of δ S,I and δ S,II . In this case there were important disparities between x TM and x Couv . 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 overestimating the contribution of the lowest layer. We therefore suggest to always attempt to fully characterize the soil isotopic profile before aggregating the isotopic information when defining the two water sources. Figure 7 gives the relative contributions from soil layers I, II, and III (upper, middle, and lower panels, respectively) to T provided by IsoSource, the multi-source mixing model of Phillips and Gregg (2003) (x IsoS , in %, displayed in the form of gray histograms), and by SIAR, the Bayesian method of Parnell et al. (2010) (x SIAR , in %, gray probability density curves). The colored vertical lines are x I_Couv , x II_Couv , and x III_Couv , the simulated S(z)dz T /( x· y) ratios from layers I, II, and III. The color-shaded areas associated with x I_Couv , x II_Couv , and x III_Couv refer to their uncertainty by accounting for the uncertainty of the input data. As for Fig. 5, δ Ti_Couv is reported above each plot along with its standard deviation. The x J _IsoS probability distribution was observed to be either narrow (e.g., for the experiment with the deep groundwater table and the dry upper layer at a low transpiration rate -DeDr_lT/layer I, Fig. 7m) or broad (e.g., for the experiment with the deep groundwater table and the wet upper layer at a high transpiration rate -DeWe_hT/layer I); i.e., the range of possible solutions for x J _IsoS was relatively small or large (10 and 100 %, respectively, for these two examples). In general, both IsoSource and SIAR results were in good agreement: the x SIAR 's most frequent value (MFV, at the peak of the density distribution curve) was in most cases either located near the median value of the x IsoS probability range (e.g., for the experiment with the shallow groundwater table and wet upper layer at a low transpiration rate -ShWe_lT/layer I, Fig. 7g) or matched exactly the x IsoS unique value (i.e., for the experiment with the deep groundwater table and dry upper layer at a low transpiration rate -DeDr_lT/layer I, Fig. 7m). In contrast, the statistical methods succeeded best in providing x estimates similar to those of the model of Couvreur et al. (2012) in the case of a shallow groundwater table and at low T only (Fig. 7a-c and g-i), i.e., when water availability was high and root compensation was low. In these cases, x I_Couv was included in the estimated x I_IsoS range and the mean absolute difference (MD) between x J _Couv and x SIAR MFV was equal to 8.6 %. This difference was greatest (129.2 %) for experiment DeDr_lT, when hydraulic redistribution by the roots was simulated by the analytical model ( Fig. 7m-o).
When considering eight soil layers instead of three, uncertainties in the assessment of the relative contributions to T as determined by IsoSource increased. The estimated probability ranges increased in most of the cases (results not shown). However, it considerably improved the results of the Bayesian model: the mean absolute difference between x J _Couv and the most frequent x SIAR value was equal to 4.7 % for the scenarios with a shallow groundwater table and low transpiration rate and equal to 52.1 % in the case of hydraulic redistribution by the roots (Table 3).
Independent of the number of defined soil layers, lowering the value of increment to 5 % in IsoSource refined the analysis where the probability distribution was already narrow (i.e., in the case of a well-identified x IsoS value, e.g., Fig. 7m), while it produced distributions that were flatter and contained fewer 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 combination and vice versa (results not shown). An increase or decrease of a factor 2 in the number of runs as well as the number of runs to be discarded from the analysis had only a marginal impact on the density distribution curves obtained with the SIAR model in the case of three or four soil layers.  (Phillips and Gregg, 2003) (displayed in the form of gray histograms). Density distribution functions following the SIAR Bayesian model (Parnell et al., 2010) (gray lines). "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. Suffices "lT" and "hT" refer to "low" and "high" transpiration rate simulations. The colored vertical lines give x I_Couv , x II_Couv , and x III_Couv , the relative contributions to transpiration from layers I, II, and III simulated by the analytical model of Couvreur et al. (2012). The color-shaded areas associated with x I_Couv , x II_Couv , and x III_Couv refer to their uncertainty associated with input data uncertainty.
The modeling 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 single or multiple "root water uptake depths", but rather a continuous RWU profile or statistical solutions of contributions to transpiration (provided by IsoSource and SIAR). Significant changes in δ Ti do not necessarily mean important changes in the depth of RWU, but rather a slight (but nevertheless significant) modification of the RWU profile. The authors believe that the relatively novel statistical tools presented in this review should therefore be preferred over the graphical inference method, especially since the two former are available as user-friendly programs and packages and do not require significant computing time, and therefore can be run locally on a personal computer. As highlighted in this series of virtual experiments, the Bayesian method showed for the case of two and three soil layers much more convincing results than the method of Phillips and Gregg (2003). The Bayesian method was particularly efficient in the case of eight soil layers, illustrating the interest of reaching the best vertical resolution and maximizing the number of identified potential sources (Table 3). Note that no prior information on the relative contributions to T from the different soil layers was used when running the SIAR Table 3. Most frequent value (mfv) and range of the density distribution curve of the relative contribution to transpiration across eight defined soil layers as determined by the Bayesian method of Parnell et al. (2010) (x SIAR , %) and the mean relative contribution (with standard deviation) provided by the analytical model of Couvreur et al. (2012) (x Couv , %). Profiles of relative contribution were computed for eight soil-plant virtual experiments differing in the depth of the groundwater table (shallow -Sh/deep -De), the soil surface water status (dry -Dr/wet -We), and the plant transpiration rate (low -lT/high -hT).

Soil layer (m)
Shallow groundwater table (Sh) Dry surface conditions (ShDr) Wet surface conditions (ShWe)  (2) 18(0-46) 23 (2) 16(0-48) 21 (2) 16(0-53) 23(2) 1.5-2 17(0-59) 14 (2) 17(0-47) 12 (2) 15(0-52) 11 (2) 16(0-51) 12 (2) Soil layer (m) Deep groundwater table (De) Dry surface conditions (DeDr) Wet surface conditions (DeWe)  (2) 18(0-53) 9(2) 16(0-46) 12 (2) program; i.e., the authors opted for flat priors. This can be changed by the user, based on additional collected data such as, for instance, information on root architecture and function across the soil profile or information on soil hydraulic properties and water status. One can show from this inter-comparison of methods that labeling 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 predominantly. However, this was never the case when considering the results of the analytical model. A dual isotope ( 18 O or 2 H) labeling 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 more accurate estimates.

Offline destructive versus online nondestructive isotopic measurements in plant and soil waters
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 resolution should, when possible, match the exponential decrease in isotopic composition (Wang et al., 2010), and it should capture sudden variations with time at the soil surface due to precipitation, i.e., be maximal at the surface and minimal deeper in the soil profile where isotopic dynamics are less pronounced. Not doing this can lead to a situation where source partitioning is not feasible from isotopic measurements. Under field conditions (i.e., ∼ 95 % of the studies reviewed in this work, summarized in Table 1) soil material is generally not a limit-Biogeosciences, 14, 2199-2224, 2017 www.biogeosciences.net/14/2199/2017/ ing factor, and thus can be sampled twice or thrice to average out or characterize lateral heterogeneity without significant disturbance of the soil. 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). Accuracy of this extraction method was shown to be maximal at higher water content and for sandy soils and lower for soils with high clay content (e.g., Koeniger et al., 2011;West et al., 2006). 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).
Certainly one of the main limitations of all isotopic approaches for quantifying RWU is the destructive character of isotopic sampling (see Sect. 3.1) and associated offline analyses (Sect. 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 measuring frequency of other soil and plant state variables, e.g., soil water content and potential, and leaf water potential (Sect. 3.2.2). In addition, one may question the representativeness of plant samples, in which tissues and thus water with very different water residence time is mixed.
Recently developed methods take advantage of laser-based spectroscopy which allows in situ, online, and continuous isotopic measurements in the gas phase at high frequency. These methods rely on coupling a laser spectrometer (e.g., wavelength-scanned cavity ring-down spectroscopy -WS-CRDS, Picarro Inc., Santa Clara, CA, USA; cavity-ringdown laser absorption spectroscopy -CRLAS; and off-axis integrated cavity output spectroscopy -ICOS, Los Gatos Research, Los Gatos, USA) with specific soil gas sampling probes consisting of gas-permeable microporous polypropylene membranes or tubing. These membranes or tubing exhibit strong hydrophobic properties, while their microporous structures allow the intrusion and collection of soil water vapor. Several authors Gangi et al., 2015;Herbstritt et al., 2012;Oerter et al., 2016;Rothfuss et al., 2013;Sprenger et al., 2015;Volkmann and Weiler, 2014) were able to 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 phases in the soil. In contrast to "traditional" isotopic methods, these novel isotopic monitoring methods also have the distinct advantage of determining soil liquid water isotopic composition at very low water content, since water vapor, in contrast to soil liquid water, is not limited 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 the laser spec-trometers could be, as pointed out by Gralher et al. (2016) for the specific case of a Picarro WS-CRDS, significantly sensitive to the carrier gas used. In their opinion papers, Mc-Donnell (2014) and Orlowski et al. (2016a) urged a comparison between methods, which was addressed by Orlowski et al. (2016b) and Pratt et al. (2016) (for vapor measurements only).
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 tested for other (nonwoody) plant species. This will imply the major 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 (Farquhar et al., 2007).

Call for a coupled experiment-modeling approach
for determination of plant water sources on the basis of isotopic data In order to fully benefit from the potential of water stable isotopic analysis as a tool for partitioning transpiration flux, the authors call for the development of approaches making use of physically based models for RWU and isotopic fractionation to analyze experimental data, especially since several soilvegetation-atmosphere transfer (SVAT) models are available that can simulate flow of isotopologues in the soil and the plant (i.e., SiSPAT-Isotope, Braud et al., 2005; Soil-litter-iso, Haverd and Cuntz, 2010;R-SWMS, Meunier et al., 2017;TOUGHREACT, Singleton et al., 2004;HYDRUS, Sutanto et al., 2012).
To the authors' knowledge, there are only a few studies which attempted to do so. Rothfuss et al. (2012) ran an experiment under controlled laboratory conditions where they measured on four dates (corresponding to four different stages of vegetation and therefore root development) soil water potential and isotopic composition profiles, and root length density distribution profiles. In their experiment, the isotopic composition of transpiration was known. The authors used a global optimization algorithm to obtain the set of parameters of SiSPAT-Isotope that best reflected the experimental dataset. Distributions of RWU could be determined on these four dates. Also, in the study of Mazzacavallo and Kulmatiski (2015), the RWU model of HYDRUS could be parameterized during a labeling (heavy water 2 H 2 O) pulse experiment on the basis of measurements of xylem water hydrogen and oxygen isotopic compositions. This provided insights into the existence of niche complementarity between tree (mopane) and grass species. Note, however, that this HY-DRUS version did not incorporate isotopic transport through the soil and the roots.
Another example is the work of Ogle et al. (2004), who were able to reconstruct active root area and RWU profiles from isotopic measurements using the 1-D analytical macroscopic model of Campbell (1991) in a Bayesian framework (root area profile and deconvolution algorithm -RAPID). By assuming normal a priori distributions for the xylem water oxygen and hydrogen isotopic compositions and considering prior knowledge on RWU distribution (i.e., synthetic information based on measurements of other studies), Ogle et al. (2004) obtained a posteriori distributions of x of a desert shrub (Larrea tridentate).
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 uncertainty, on the condition that all required (isotopic) data are 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 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 field -but 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 the root function (active versus non-active roots in the soil profile) in water uptake and, thus, quantify the disconnection between measured RLD and the prognostic variable SSF (see Appendix B1). To do this, controlled conditions in state-ofthe-art climatic chambers are ideal, as they allow for a reduction of 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 also impose lower boundaries for the soil compartment (e.g., drainage and capillary rise dynamics) and provide the means to close the hydrological balance are 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 the mimicking of natural conditions.

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 RWU 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 a powerful and valuable tool for the assessment of plant water sources. In an inverse modeling framework, water isotopic analysis of plant tissues and soil also allows for obtaining of species-specific parametrization of physically based analytical and numerical RWU models. They provide a unique way to tackle the difficulty of disentangling actual RWU profiles with root traits and characteristics. Yet our literature review revealed that isotopic data have been up to now mainly used to assess water sources under natural ecosystems using statistical approaches. Only 4 % of current scientific publications demonstrate the use of a physically based model for analyzing isotopic data.
Three methods (representing 90 % of the studies) have been used for partitioning water sources: the "direct inference" method, the two-end-member mixing model and two examples of multi-source mixing models. We performed a comparison between these methods. We were able to quantify the impact of the definition of the plant water sources (i.e., whether they are spatially disjoint or not and whether their isotopic composition values are significantly different or not) on the outcome of the two-end-member mixing model. We highlighted the importance of systematically reporting uncertainties along with estimates of contribution to T of given plant water sources. 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 good performance of the Bayesian approach of Parnell et al. (2010), which uses a more rigorous statistical framework, if the number of considered water sources matches the number of isotopic measurements in the soil profile. However, in contrast to the analytical model, none of the graphical and statistical methods was able to locate and quantify hydraulic redistribution of water.
Finally, the authors call for a generalization of coupled approaches relying on the confrontation between labeling experiments under controlled conditions and three-dimensional RWU numerical modeling. This type of approach could be used in agronomy to quantify RWU as a function of plant genotype and soil structure. It also has great potential for quantifying RWU in seminatural and natural ecosystems for understanding the mechanisms underlying the vegetation feedbacks to the atmosphere in the contexts of land cover and climate changes.
Data availability. Data used by the authors are available in Table 2 x j , x Couv , x J _Couv , x J _IsoS , x J _SIAR Relative contribution to transpiration, source j relative contribution to transpiration, continuous and integrated (layer J ) relative contributions to transpiration as simulated by the analytical model of Couvreur et al. (2012), integrated (layer J ) relative contributions to transpiration as determined by the IsoSource and SIAR models of Phillips and Gregg (2003) and Parnell et al. (2010). Relative contribution to transpiration under conditions of uniform soil water potential -7, 8b, 9, 9 s z, z j , z j +1 , z j , z max , z RWU 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" L 4b, 5, B2-B4 m/p δ, δ 2 H, δ 18 O, δ source , δ surf , δ S , δ S,j , δ S,J , δ A , δ B , δ Ti , δ Ti,m , δ Ti_Couv , δ E Water stable isotopic composition, water hydrogen and oxygen stable isotopic compositions, source, soil surface, soil water, 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) -(expressed in ‰) 3-5, 7-9, 11 m/s ε j Residual error term -(expressed in ‰) 9 s θ Soil volumetric water content L 3 L −3 5, 11 m ρ Volumetric mass of water M L −3 3 m σ x , σ δ A , σ δ B , σ δ Ti , σ δ Ti_Couv , σ x Standard errors associated with the measurements of x, δ A , δ B , and δ Ti and estimated uncertainty of δ Ti_Couv as simulated by the analytical model of Couvreur et al. (2012), error associated with the estimation of the relative contribution to T of water source A in the case of two distinct sources -(expressed in ‰) 8a, 8b s τ Isotopic tolerance -(expressed in ‰) 10 p