Reviews and syntheses: Gaining insights into evapotranspiration partitioning with novel isotopic monitoring methods

Disentangling ecosystem evapotranspiration (ET) into evaporation (E) and transpiration (T) is of high relevance for a wide range of applications, from land surface modelling to policymaking. Identifying and analysing the determinants of the ratio of T to ET ($T/ET$) for various land covers and uses, especially in view of climate change with an increased frequency of extreme events (e.g. heatwaves and floods), is prerequisite for forecasting the hydroclimate of the future and tackling present issues, such as agricultural and irrigation practices.

One partitioning method consists of determining the water stable isotopic compositions of ET, E, and T (δ_ET, δ_E, and δ_E, respectively) from the water retrieved from the atmosphere, the soil, and the plant vascular tissues. The present work emphasizes the challenges this particular method faces (e.g. the spatial and temporal representativeness of the $T/ET$ estimates, the limitations of the models used, and the sensitivities to their driving parameters) and the progress that needs to be made in light of the recent methodological developments. As our review is intended for a broader audience beyond the isotopic ecohydrological and micrometeorological communities, it also attempts to provide a thorough review of the ensemble of techniques used for determining δ_ET, δ_E, and δ_E and solving the partitioning equation for $T/ET$.

From the current state of research, we conclude that the most promising way forward to ET partitioning and capturing the subdaily dynamics of $T/ET$ is by making use of non-destructive online monitoring techniques of the stable isotopic composition of soil and xylem water. Effort should continue towards the application of the eddy covariance technique for high-frequency determination of δ_ET at the field scale as well as the concomitant determination of δ_ET, δ_E, and δ_E at high vertical resolution with field-deployable lift systems.

1 Introduction

A pivotal parameter in landscape hydrology and ecology is the transpiration (T) to evapotranspiration (ET) ratio ($T/ET$) (see the reviews of Kool et al., 2014; Anderson et al., 2017; Stoy et al., 2019). Isolating the T flux in ET is of utmost importance for a wide range of applications because of its link to plant water uptake, for e.g. optimizing irrigation practices (Skaggs et al., 2010), tackling ecological questions in water-limited ecosystems (Rothfuss and Javaux, 2017), or a better representation of the relations between the carbon and water cycles in climate models (Humphrey et al., 2018; Ito and Inatomi, 2012). At the global scale, the uncertainty of the $T/ET$ estimate remains high; it has been estimated to range from 13 % to 90 %, depending on the source and type of data (e.g. satellite- or isotopic-based) and method (modelling or data reanalysis) (Lawrence et al., 2007; Alton et al., 2009; Jasechko et al., 2013; Wang et al., 2014; Wei et al., 2017). Ultimately, this conditions the ability of land surface models to provide sensitivity of the overall ET flux to changes in precipitation and land cover (Wang and Dickinson, 2012).

Spatial and temporal variability add even more uncertainty to our knowledge about $T/ET$ at the local scale, which is a prerequisite for a meaningful use of such estimates for any of the practical and scientific questions mentioned above. Partitioning ET into the raw components E and T at the field and subdaily spatiotemporal scales is generally performed by an ensemble of partitioning methods, which can be divided into instrumental approaches and correlation-based modelling approaches (Scanlon and Kustas, 2010). The former approach includes e.g. the eddy covariance (EC) technique (Baldocchi, 2014; Reichstein et al., 2005), soil flux chamber measurements (Raz-Yaseef et al., 2010; Lu et al., 2017), micro-lysimeter measurements (Kelliher et al., 1992), or atmospheric-profile measurements (Ney and Graf, 2018).

Another instrumental method to partition ET is based on the analysis of its hydrogen or oxygen isotopic composition, i.e. the water vapour atom ratio in rare (²H or ¹⁸O) and abundant (¹H or ¹⁶O) stable isotopes and expressed on the international “delta” (δ) scale (Dubbert and Werner, 2019). The method utilizes the natural discrepancies in isotopic composition of the ecosystem evaporation (δ_E) and transpiration (δ_E) fluxes. The difference of δ_E − δ_E originates primarily from thermodynamic and kinetic fractionation during phase change and transport processes undergone by water evaporating from soil on the one hand and water extracted by a root system and transpired by the canopy on the other hand. The observed discrimination against stable isotopologues along the soil–plant–atmosphere water path can be conceptualized twofold, i.e. phase-change- and diffusion-driven, and quantified by the so-called equilibrium and kinetic fractionations, respectively, for which we will later review the physically based expressions. The term (δ_E − δ_E) is also determined by (see Fig. 1)

• i.

  the difference in boundary conditions acting on E and T, i.e. the δ value of soil water at the evaporating front (EF) and that of the leaf water at the transpiration site and of the atmospheric water vapour;

• ii.

  the prevalence (or non-prevalence) of isotopic steady state (ISS) for transpiration, i.e. whether δ_E is independent of time (Farquhar and Cernusak, 2005; Dubbert et al., 2014a) (see Sect. 3 for a detailed description of ISS; note that the ISS assumption is generally not made for evaporation flux, but see Rothfuss et al., 2010, for an exception).

The spatiotemporal variabilities of these factors and the complexity of their interactions may result in significant heterogeneous distributions of both δ_E and δ_E in the field (Fig. 1). Importantly and as reflected by the reviewed isotopic literature (see Sect. 2), E in this context does not include canopy interception and dew evaporation, which are known to be associated with isotopic effects (Allen et al., 2017; Zheng et al., 2019). Theses fluxes can be of significant magnitude, depending on the scale of interest (Good et al., 2015; Allen et al., 2017). The $T/ET$ fraction is obtained by inverting the isotopic mass balance equation ${\mathit{\delta }}_{ET}=(\mathrm{1}-T/ET){\mathit{\delta }}_E+(T/ET){\mathit{\delta }}_{\mathrm{T}}$ (Yakir and Sternberg, 2000):

$$\begin{array}{}\text{(1)}& T/ET={\displaystyle \frac{{\mathit{\delta }}_{ET}-{\mathit{\delta }}_E}{{\mathit{\delta }}_{\mathrm{T}}-{\mathit{\delta }}_E}}.\end{array}$$

Equation (1) highlights how the isotopic partitioning methodology differs from other instrumental approaches, such as those based on a combination of different techniques (e.g. lysimeter and EC measurements): it solely relies on measurements and/or analytical modelling of the stable isotopic compositions of the components ET, E, and T. Behind this apparent simplicity and the problem of (e.g. spatial) representativeness highlighted in Fig. 1 put aside, the isotopic partitioning methodology is limited in its application in different ways, such as the inability – until recently – to provide continuous (i.e. non-destructive) δ_E, δ_E, and δ_ET assessments. Part of these limitations were overcome with the availability of field-deployable laser-based spectrometers. These instruments allow for long-term monitoring of soil water vapour and plant transpiration isotopic compositions when combined with gas-permeable membrane or tubing technology (Beyer et al., 2020).

A variety of different methods exist to measure or estimate δ_E, δ_E, and δ_ET. The central aim of this study is to identify from the literature the challenges the ensemble of isotopic methods currently face and how they should progress in the future (Sect. 3). Particularly, the abovementioned emerging monitoring methods are reviewed for the specific purpose of ET partitioning. As such, our work differs from those of Wang and Yakir (2000), Yakir and Sternberg (2000), Xiao et al. (2018), and Sun et al. (2019). Note also that this study will not focus on differences in $T/ET$ as estimated by the abovementioned traditional methods on the one hand and by the isotopic methods on the other; this has been extensively reported by e.g. Sutanto et al. (2014). In addition and for non-specialist readers, we thoroughly review the underlying concepts and techniques involved in the determination of δ_E, δ_E, and δ_ET. In order to highlight the important progresses made over the past 30 years, we also give a literature overview (Sect. 2). Finally, Sect. 5 presents a summary of our own suggestion for improvement as well as of the possible ways forward for the isotopic partitioning community.

Figure 2 Graphical summary of the reviewed literature. (a) Evolution of the number of citations per year (blue bars) and cumulative number of publications (1989–2020, red line). (b) Temporal resolution vs. extent of the estimate of the transpiration-to-evapotranspiration ratio ($T/ET$). Numbers above or below the histograms refer to the number of studies working at a given temporal resolution. (c, d) Listing of the different plant cover and climate types with proportions (white label) expressed in percentage and (g) map locating each study (with reference number 1–39 with a white label; see Table A1). (e, f) Proportions of field vs. modelling studies and prevailing experimental conditions (as natural precipitation or irrigation or else as labelling studies). (h–j) Listing and proportions of methods for determination of δ_ET, δ_E, and δ_E.

2 Literature overview

A total of 39 studies were found by entering the search terms ((“evapotranspiration” or “transpiration” or “evaporation”) and partition* and isotop*) into the ISI (Institute for Scientific Information) Web of Science search engine (http://www.webofknowledge.com/, last access: 15 February 2021). The reader will find a graphical summary in Fig. 2 as well as a detailed description for each of the entries in Table A1 of Appendix A. On average, approximately 1.3 (2.4) partitioning studies were published each year over the period 1989–2007 (2008–2020) with an average annual citation rate of 12 (143) (Fig. 2a).

To the authors' knowledge, the first scientific article reporting on the possibility to partition ET on basis of the differences in isotopic composition of ecosystem ET, soil evaporation, and plant transpiration was that of Bariac et al. (1987). An attempt to use this possibility was made in the study of Walker and Brunel (1990) (Table A1) but remained, according to the authors, not conclusive; 10 years later, Jean-Pierre Brunel and his colleagues could provide the first water-stable-isotope-derived estimation of the relative importance to ET of the transpiration of the tropical and water-stress-resistant plant Guiera senegalensis (Brunel et al., 1997), which was noticeably low (approx. 20 %). In the meantime, Moreira et al. (1997) applied the so-called “Keeling plot” technique (Keeling, 1958) (see Sect. 3.1) for determination of the isotopic composition of ET for the specific purpose of partitioning. The isotopic compositions of soil E and plant T at two sites (one pasture and one forest) in the Amazon basin were inferred by using the atmospheric part of the Craig and Gordon (1965) model (see Sect. 3.2) and by assuming steady-state transpiration flux (see Sect. 3.3), respectively. The authors could provide evidence of the strong prevalence of T in the ET budget. In a hybrid work coupling a review of the state of the art with field measurements, Wang and Yakir (2000) concluded on the predominance of T flux in a wheat field located in the Negev region, Israel (i.e. $T/ET$ > 96.5 %).

As represented in Fig. 1, partitioning ET may be significantly complicated in cases of mixed vegetation covers. A few studies focused on estimating the ratio of the vegetation type or strata-specific transpiration to evapotranspiration. Yepez et al. (2003) applied the Keeling plot technique specifically to two distinct ecosystem layers of a savanna woodland in southern Arizona, USA, i.e. the understorey dominated by the Sporobolus wrightii C₄ grass and the canopy populated by the mesquite tree Prosopis velutina. By doing this, they could capture the isotopic composition of ET representative of each of the two ecosystem layers. In order to partition ET, the authors computed the isotopic composition of the whole ecosystem T as a composite function of the isotopic compositions of grass and tree T fluxes. Finally, it was determined that grass and tree T amounted to 15 % and 75 % of total ET. Xu et al. (2008) investigated the discrepancies between $T/ET$ assessments from either δ²H or δ¹⁸O data collected in a subalpine shrubland (Balang Mountain, China). They could differentiate between tree (Quercus aquifolioides) and understorey (e.g. Cystopteris montana) contributions to ET by using the multi-source mixing model IsoSource (Phillips and Gregg, 2003). In an open cork-oak (Quercus suber L.) savanna, Dubbert et al. (2014b) investigated the impact of the understorey vegetation (annual grass and forbs) on the total ecosystem water budget. They could discriminate between T of trees and grass and highlighted the stability of the former throughout the year and the strong decrease of the latter during the summer. Piayda et al. (2017) differentiated between open and shaded portions of the same experimental site and found $T/ET$ ranging from 9 % to 59 % and from 17 % to 66 %, respectively. Zhang et al. (2018) investigated a marsh wetland in China and found out that the two dominant species (Scirpus triqueter and the invasive Phragmites australis) contributed equally (20 %) to ET flux.

A number of authors either investigated the impact of irrigation on the partitioning of ET or relied on irrigation pulses, i.e. applied volumes of isotopically enriched or depleted water (with respect to local irrigation water) to the soil. By doing this, they could reduce the uncertainty of the $T/ET$ estimates by artificially enhancing the difference between δ_E and δ_E. In a study conducted in a semi-arid environmental setting (Marrakesh, Morocco), Williams et al. (2004) observed that irrigation enhanced soil E of an olive orchard (Olea europaea L.). Midday average $T/ET$ decreased from approx. 100 % (determined prior irrigation) to 69 %–86 % (computed over the 5 d period after irrigation). Yepez et al. (2005) used large gas exchange chambers either positioned on bare-soil plots or sparsely vegetated areas of a semi-arid grassland in Arizona, USA. They determined with the Keeling plot technique the isotopic composition of E and ET following an irrigation pulse. This is, to the authors' knowledge, the first use of a closed chamber system in the context of ET partitioning where T is the single source of the change in air moisture concentration. In contrast to the previous partitioning studies, Yepez et al. (2005) determined the isotopic composition of the non-steady-state (NSS) T flux on the basis of plant physiological and micrometeorological measurements using the formulation of Farquhar and Cernusak (2005) (see also for later examples Sun et al., 2014; Hu et al., 2014). The authors finally calculated $T/ET$ values ranging between 35 % and 43 % the first 3 d after irrigation, and these decreased to 22 % after 1 week. Aouade et al. (2016) found a decreasing diurnal (i.e. morning vs. afternoon) amplitude of $T/ET$ in a winter wheat field in Morocco under wet conditions after flood irrigation (soil water content value of approx. 0.35 m³ m⁻³) and the opposite under dry conditions (soil water content value of approx. 0.15 m³ m⁻³). Aouade et al. (2020) compared the $T/ET$ results for dry conditions of Aouade et al. (2016) to independent assessments using the Interaction between Soil, Biosphere, and Atmosphere (ISBA) model (Masson et al., 2013) and found that they were within the same range (73 %–89 %). In another study, Good et al. (2014) found on average a value of 30 (± 5) % for $T/ET$ in a grassland site during the first 15 d following a 30 mm isotopically enriched irrigation event. Finally, Lu et al. (2017) focused on the efficiency of irrigation strategies in southern California (USA). They documented that the investigated field of Sorghum bicolor was responsible for 46 % of water consumption following the irrigation event.

Hsieh et al. (1998), Ferretti et al. (2003), Wenninger et al. (2010), and Sutanto et al. (2012) obtained $T/ET$ values by the closing of a common water isotope mass balance equation. For this, the authors made a series of simplifying hypotheses: atmospheric water vapour is in thermodynamic equilibrium with soil water, and the isotopic composition of T is the amount-weighted average of the isotopic compositions of precipitation and soil water. Ferretti et al. (2003) obtained $T/ET$ values ranging between 10 % and 60 %, depending on the growing season, in a semi-arid grass steppe, while Hsieh et al. (1998) estimated $T/ET$ to span from 14 % to 71 % as annual rainfall increased along two sampling transects in Hawaii. Wenninger et al. (2010) and Sutanto et al. (2012) applied the isotope mass balance equation in similar semi-controlled experimental setups equipped with soil liquid water (Rhizon) samplers. In their framework, the destructive sampling of soil to retrieve the isotopic composition of soil E was not needed, while a number of simplified hypotheses had to be made regarding T. Wenninger et al. (2010) simulated a $T/ET$ value of 70 % for teff (Eragrostis tef) during the course of their experiment. Sutanto et al. (2012) found a comparable value for a grass cover ($T/ET$ = 87 %). In both of these studies, the isotopic partitioning results were confronted with additional (e.g. micrometeorological) measurements and independent models such as HYDRUS-1D (Simunek and van Genuchten, 2008).

Isotope-enabled, physically based, and numerically solved soil–vegetation–atmosphere transfer (SVAT) models were also tested against $T/ET$ data collected in both laboratory and field setups. In the study of Rothfuss et al. (2012), $T/ET$ of a 0.2 m² surface area monolith was simulated with the SiSPAT-Isotope model (Simple Soil–Plant–Atmosphere Transfer; Braud et al., 2005) at five selected dates under strictly controlled conditions in a climate chamber along the development of a tall fescue cover (Festuca arundinacea). $T/ET$ was determined to increase from 6 % (16 d after sowing) to 95 % (43 d after sowing); 1 year earlier, Haverd et al. (2011) used another isotopically SVAT model, Soil–Litter–Iso (Haverd and Cuntz, 2010), using data from a field experiment (eucalyptus forest, southeastern Australia) in a similar framework, i.e. by running a multi-objective calibration to estimate a given set of model parameters. However, in contrast to Rothfuss et al. (2012), they could show that the added information provided by the isotopic data (δ²H) was not effective in better constraining the model for determination of $T/ET$ (in their case equal to 85 ± 2 %). Another simulation study was published by Wang et al. (2015), where a physically based model solving the energy and water balance in the soil–plant–atmosphere continuum (Wang and Yamanaka, 2014) was coupled to an isotopic module accounting for fractionation processes during E and T. Wang et al. (2015) simulated $T/ET$ of a temperate grassland to spread over a wide range of values (i.e. 2 %–99 %) during the course of a 190 d long experiment. Wei et al. (2018) used a similar modelling framework as in Wang et al. (2015) and found that the 3-month ET-weighted $T/ET$ values were equal to 74, 93, and 81 % for three different crops, i.e. rice, corn, and wheat, respectively, grown in temperate (rice, Japan) and semi-arid monsoonal (corn and wheat, China) environmental conditions.

Wang et al. (2010, 2013) published the first ET partitioning studies where water vapour hydrogen and oxygen isotopic compositions were measured online with an infrared laser spectrometer. Using closed gas exchange chambers, they determined by mass balance the isotopic compositions of E, T, and ET in a non-destructive way. This allowed the authors not to rely on either (i) making the assumption of T at ISS for partitioning ET fluxes or (ii) modelling the isotopic composition of T at NSS (see Sect. 3.3). Wang et al. (2010) calculated $T/ET$ values for the mesquite tree (Prosopis chilensis) grown under controlled conditions (Biosphere 2 facility, Arizona, USA; see for details Barron-Gafford et al., 2007), ranging from 61 % to 83 % at 25 % and 100 % woody cover, respectively. Wang et al. (2013) compared $T/ET$ values (65 %–77 % vs. 83 %–86 %) computed from control vs. warming plots, taking advantage of a long-term grassland multiple-factor climate control experiment in Oklahoma, USA.

The quantification of the overall uncertainty associated with isotope-derived $T/ET$ estimates has been the focus of several studies. Other studies focused on the sensitivity of $T/ET$ to different environmental (e.g. isotopic) factors. Good et al. (2014) studied for instance the uncertainty of the $T/ET$ values obtained at their grassland site as a function of the uncertainty linked with the estimate of δ_ET obtained with the Keeling plot technique (according to Good et al., 2012). Bijoor et al. (2011) highlighted the high uncertainty of their $T/ET$ isotope estimates (i.e. standard error value > 37 %). Xu et al. (2008) and Yepez et al. (2005) calculated the uncertainties linked with determination of $T/ET$ with the IsoError software (Phillips and Gregg, 2001). Dubbert et al. (2013) quantified the sensitivity of the partitioning of ET to a number of factors (e.g. value of the kinetic fractionation factor and assumption of steady-state T) during a field experiment in central Portugal. They also compared direct measurements of the isotopic composition of E (with gas exchange chambers coupled to a laser spectrometer) to simulations with the Craig and Gordon (1965) model. Similar to Rothfuss et al. (2010, 2012), the authors underlined the need to complement isotopic measurement with micrometeorological and physiological observations. Hu et al. (2014) determined a mean $T/ET$ value of 83 % in a semi-arid shrubland in China dominated by Stipa krylovii and Artemisia frigida. They tested for the first time the so-called flux gradient approach (Lee et al., 2007; see Sect. 3.1) for determination of δ_ET. The authors argued that, in their case, the uncertainty of the δ_ET estimates had the strongest effect on $T/ET$ uncertainty. Also Wei et al. (2015) found that the greatest source of uncertainty of $T/ET$ of a rice paddy field was linked to the determination of δ_ET, this time using the Keeling plot technique. They could further express $T/ET$ as an exponential function of leaf area index (LAI) (i.e. $T/ET$ = 67 ⋅ LAI^0.25, expressed in %) at the seasonal scale. Wu et al. (2017) found slightly different parameters of the same LAI model (71 ⋅ LAI^0.14) for a maize crop grown under semi-arid conditions (Gansu Province, China).

Among the studies listed in Table A1, a few complemented their isotopic methods with traditional instrumental approaches, such as EC, soil flux chambers, and lysimeters, and investigated the goodness of fit between the isotopic and non-isotopic $T/ET$ values. Sutanto et al. (2014) reported from the literature generally higher isotope-derived $T/ET$ (> 70 %) values than those of the traditional approaches for comparable land cover types. However, at experimental sites combining both type of measurements, Sutanto et al. (2014) underlined a fair agreement between both approaches. Bijoor et al. (2011) investigated the partitioning of ET in a freshwater marsh dominated by Typha latifolia in California, USA. They found a good agreement between $T/ET$ values estimated on the one hand from isotopic analysis and from micrometeorological (e.g. EC) measurements on the other. Berkelhammer et al. (2016) compared the outcome of the isotopic partitioning method with EC-derived $T/ET$ values. They underlined the goodness of fit of the two methods as well as the stability of $T/ET$ as a function of LAI over multiannual timescales. Wen et al. (2016) investigated the contribution of spring maize T to ET in an arid artificial-oasis part of the Heihe River catchment (China) and reported it to be quite constant (mean $T/ET$ value of 87 ± 5.2 %). Collected data were further used by Zhou et al. (2018) and Xiong et al. (2019). Zhou et al. (2018) showed similarities between results of the isotopic partitioning method and a coupled approach of EC and lysimeter data. They underlined, however, that both methods simulate higher $T/ET$ values, with poor temporal dynamics not reflecting those of the leaf area index, than their benchmark method, i.e. based on the incorporation of the vapour pressure deficit into the expression of the water use efficiency concept. Xiong et al. (2019) observed a good match between $T/ET$ daily values (54 %–97 %, with a mean value of 85 %) as obtained with their isotope method and with a net radiation and temperature-dependent model coupled to imaging radiometry. Quade et al. (2019) cross-compared the $T/ET$ values based on either water δ²H or δ¹⁸O data at selected dates along the development of a sugar beet (Beta vulgaris) crop with different methods including the combination of EC and lysimeter flux data.

Until now, only two studies have made use of gas-permeable membranes for online and non-destructive determination of δ_E and determination of $T/ET$ values. Gaj et al. (2016) fitted a one-dimensional analytical solution of the water isotopic composition in the soil profile to their data to retrieve $T/ET$ values in the semi-arid Cuvelai-Etosha Basin, Namibia. Quade et al. (2019) compared $T/ET$ results obtained on the basis of the non-destructive method of Rothfuss et al. (2013) with those of traditional destructive soil sampling. They found significant differences in $T/ET$ between the different methods on 4 d at different stages of the sugar beet canopy development (0.7 < LAI < 6.7).

In a review on the use of water stable isotope analysis for determination of plant root water uptake dynamics (Rothfuss and Javaux, 2017), the authors underlined the need for field studies in croplands. This is not the conclusion of the present literature overview, as the three main land surface types, i.e. cropland, forests, and grassland (in monoculture or mixed culture), are rather equally represented with a relative proportion of 33 %, 32 %, and 41 %, respectively (Fig. 2b). More than one-third of the scientific publications analysed in the present review (i.e. 38 %) applied the isotopic methodology in semi-arid or desert regions (Fig. 2c). Nevertheless, a wide range of climate types (e.g. subtropical humid, Mediterranean, or subarctic, Fig. 2c) as well as regions (e.g. North America, sub-Saharan Africa, or eastern Asia, Fig. 2d) is investigated as well. Of the 39 reviewed studies 30 were conducted in the field, and only 8 (21 %) used a physically based numerical model to simulate $T/ET$ on the basis of the collected isotopic (and water status) data (Fig. 2e). Furthermore, 95 % of the field studies were conducted at natural isotopic abundance, either under a normal precipitation regime (85 %) or in the framework of an irrigation experiment (10 %). The remaining 5 % of studies (Yepez et al., 2005; Good et al., 2014) applied a labelling pulse of ²H-enriched water to the soil for better discrimination between the three terms of the mixing equation (Eq. 1).

Figure 3 Summary of the different approaches (mass balance, physically, and statistically based) methods for determination of δ_ET with the relevant variables and fluxes for each case.

There is naturally a strong link between the temporal resolution in $T/ET$ estimates and the temporal extent of the $T/ET$ time series (Fig. 2b). The vast majority of the studies (85 %) provided $T/ET$ values at hourly to subweekly resolutions over periods of time not exceeding a few months. This is partly a sign of the limitation of the isotopic methodology, which was mentioned in the Introduction, i.e. the labour-intensive and time-consuming destructive sampling of soil and plant material and the subsequent water extraction step. In two studies only (Hsieh et al., 1998; Ferretti et al., 2003), authors could calculate $T/ET$ at weekly to monthly resolutions over several years. For doing this, they made a series of abovementioned simplifying hypotheses, which allowed them, amongst other things, not to rely on sampling of plant material, thereby significantly saving extraction and analysis time. The authors of the present work note that, on the other hand, the question of spatial variability or representativeness of the $T/ET$ estimates is rarely addressed in the literature (but see Sect. 3.1 for the issue of spatial representativeness of δ_ET).

3 Methodological review

In this section, the methods used for determination of the three terms in the partitioning equation (Eq. 1), i.e. δ_ET, δ_E, and δ_E, for final computation of $T/ET$ will be covered (Sect. 3.1.1, 3.2.1, and 3.3.1, respectively), with special emphasis on challenges and new technical and methodological developments (Sect. 3.1.2, 3.2.2, and 3.3.2, respectively). Three main approaches emerge from the analysis (Fig. 3): δ_ET, δ_E, and δ_E can be either determined by

• i.

  solving the mass balances for the different water vapour isotopologues,

• ii.

  using physical models based on macroscopic analogies of Ohm's law, or

• iii.

  using a statistical framework.

Note that it is not the present work's intention to give a thorough review of the physically based and isotope-enabled soil–vegetation–atmosphere numerical models used by Haverd et al. (2011), Sutanto et al. (2012), and Rothfuss et al. (2012) for simulation of $T/ET$. For this, the readers may refer also to Haverd and Cuntz (2010) and Braud et al. (2005). Likewise, the authors choose not to describe one particular ensemble of methods in detail (used in seven different studies; see Table 2 and referred to as “water balance” in Fig. 2.) based on solving a water mass balance equation and not relying on the sampling or monitoring of plant and soil water and atmospheric water vapour.

3.1 Isotopic composition of evapotranspiration

3.1.1 Methods

The prevalent method (43 % of the reviewed studies, Fig. 2h) for determining the isotopic composition of ET is based on solving a mass balance equation (Fig. 3a–c). It was named after Charles D. Keeling who originally used it to quantify the CO₂ carbon isotopic composition in the atmosphere as a linear function of the reciprocal of the CO₂ concentration (Keeling, 1958). The so-called Keeling plot technique simply considers that the water vapour measured in some ecosystem atmosphere (of concentration C_atm and dimension of M L⁻³), e.g. within or above the canopy, originates from two sources, namely (i) the background water vapour (of concentration C_bg, M L⁻³), transported advectively and defined as not being influenced by ET flux, and (ii) evapotranspiration ET (of concentration C_ET, M L⁻³):

$$\begin{array}{}\text{(2)}& C_{\text{atm}}=C_{\text{bg}}+C_{ET}.\end{array}$$

Practically, laser-based spectrometers measure the water vapour volume mixing ratio, χ (–), the ratio of water vapour pressure and total (dry) atmospheric pressure:

$$\begin{array}{}\text{(3)}& {\mathit{\chi }}_{\text{atm}}={\mathit{\chi }}_{\text{bg}}+{\mathit{\chi }}_{ET}.\end{array}$$

Figure 4 Illustration of the Keeling (1958) plot technique for determination of the isotopic composition of the surface flux, here evapotranspiration (δ_ET). Subscript “bg” refers to the atmospheric background air, i.e. the air, which is not influenced by the surface ET flux. (a) Cases with different slopes of the regression line and implications for the nature of the surface flux: ET tends either toward transpiration (T) or evaporation (E). Illustration of the importance of the (narrow or wide) spread in values of the water vapour mixing ratio (χ_atm, ppmV) for the uncertainty of the δ_ET estimate (b, c).

A similar equation can be written for stable isotopes:

$$\begin{array}{}\text{(4)}& {\mathit{\delta }}_{\text{atm}}{\mathit{\chi }}_{\text{atm}}={\mathit{\delta }}_{\text{bg}}{\mathit{\chi }}_{\text{bg}}+{\mathit{\delta }}_{ET}{\mathit{\chi }}_{ET},\end{array}$$

with δ_atm and δ_bg being the isotopic compositions of the ambient air and background air, respectively. Combining Eqs. (3) and (4) and rearranging for δ_atm leads to the following expression (Eq. 5; see Fig. 3 for an illustration):

$$\begin{array}{}\text{(5)}& {\mathit{\delta }}_{\text{atm}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\chi }}_{\text{atm}}}}\left[{\mathit{\chi }}_{\text{bg}}\left({\mathit{\delta }}_{\text{bg}}-{\mathit{\delta }}_{ET}\right)\right]+{\mathit{\delta }}_{ET}.\end{array}$$

To the conditions that

• i.

  both χ and δ values of the background water vapour and ET remain constant during the measurement period and

• ii.

  there is no loss of water vapour from the atmosphere (e.g. during dewfall),

it is possible to determine δ_ET as the y intercept of the regression line of the relationship between δ_atm and $\mathrm{1}/{\mathit{\chi }}_{\text{atm}}$. In this framework the sign of the linear regression slope s = ${\mathit{\chi }}_{\text{bg}}\left(\mathrm{M}{\mathrm{L}}^{-\mathrm{3}}\right)({\mathit{\delta }}_{\text{bg}}-{\mathit{\delta }}_{ET})$ is therefore constrained by the difference (δ_bg − δ_ET); s is generally negative, apart from some bare-soil situations (Yakir and Sternberg, 2000) (Fig. 4a). Note that it is also possible to derive δ_ET by inverting the expression for s, although, to our knowledge, such a possibility has not yet been tested in the literature, certainly because the determination of C_bg and δ_bg is not straightforward in the field.

One important prerequisite for the applicability of the Keeling plot is a significant span in χ_atm values over the course of the measurements (Fig. 4b and c). High χ_atm values are especially needed to reduce the statistical uncertainty of δ_ET (Good et al., 2012). In the case of a single observation height (Wei et al., 2018, 2015; Good et al., 2014), the time factor is critical. χ_atm variations should not be obtained at the expense of the validity of the aforementioned core assumption (i), i.e. steady values of δ_ET and background χ and δ. Another option beside the one just described, which we could refer to as the “temporal Keeling plot”, is to drastically increase the span of χ_atm values by collecting data at different observation heights during a short period of time (∼ 1 h), which could be referred to as the “spatial Keeling plot”. From our literature compilation, the spatial Keeling plot is preferred over the temporal one (i.e. 32 vs. 7 studies).

Another technique (18 % of the reviewed studies, Fig. 2h) for determining δ_ET requires the manipulation of transparent chamber systems to enclose and tightly seal the soil and vegetation (e.g. Yepez et al., 2005; Piayda et al., 2017). Two different applications exist, both based on the mass balance approach. In the first one (referred to as “Chamber (InOut)” in Table A1 and Fig. 2h), the chamber is flushed with ambient air, and δ_ET is deduced from the difference in the water vapour mixing ratio and isotopic composition measured alternatingly at the inlet (subscript “in”) and outlet (subscript “out”) of the chamber (e.g. Wang et al., 2013; Dubbert et al., 2013):

$$\begin{array}{}\text{(6)}& {\mathit{\delta }}_{\text{out}}{\mathit{\chi }}_{\text{out}}={\mathit{\delta }}_{\text{in}}{\mathit{\chi }}_{\text{in}}+{\mathit{\delta }}_{ET}{\mathit{\chi }}_{ET}.\end{array}$$

Equation (6) is strictly valid only for conservative flow conditions. In other studies (e.g. Dubbert et al., 2014b), the change in flow rate (u  L³ T⁻¹) between the in- and outlet due to the addition of water vapour originating from the soil and/or the plant is taken into account as follows:

$$\begin{array}{}\text{(6')}& {\mathit{\delta }}_{\text{out}}{\mathit{\chi }}_{\text{out}}{u}_{\text{out}}={\mathit{\delta }}_{\text{in}}{\mathit{\chi }}_{\text{in}}{u}_{\text{in}}+{\mathit{\delta }}_{ET}{\mathit{\chi }}_{ET}(u_{\text{out}}-u_{\text{in}}).\end{array}$$

By conservation of dry airflow, i.e. $u_{\text{out}}(\mathrm{1}-{\mathit{\chi }}_{\text{out}})=u_{\text{in}}(\mathrm{1}-{\mathit{\chi }}_{\text{in}})$ (Simonin et al., 2013), Eq. (6') becomes

$$\begin{array}{}\text{(7)}& {\mathit{\delta }}_{ET}={\displaystyle \frac{{\mathit{\chi }}_{\text{out}}{\mathit{\delta }}_{\text{out}}-{\mathit{\chi }}_{\text{in}}{\mathit{\delta }}_{\text{in}}}{{\mathit{\chi }}_{\text{out}}-{\mathit{\chi }}_{\text{in}}}}-{\displaystyle \frac{{\mathit{\chi }}_{\text{in}}{\mathit{\chi }}_{\text{out}}({\mathit{\delta }}_{\text{out}}-{\mathit{\delta }}_{\text{in}})}{{\mathit{\chi }}_{\text{out}}-{\mathit{\chi }}_{\text{in}}}}.\end{array}$$

The second term on the right-hand side of Eq. (7) therefore accounts for the increase of flow rate due to ET in the chamber. An alternative consists of flushing the chamber with dry air instead of ambient air so that the isotopic composition of the outlet water vapour directly reflects that of ET. In the second application (named “Chamber (Keeling plot)” in Fig. 2h), the chamber is flushed in a closed loop with ambient air, and δ_ET is obtained by linear regression of the isotopic composition of the chamber air vs. the inverse of the water vapour mixing ratio using the Keeling (1958) plot technique.

In 10 % of the referenced studies (Wen et al., 2016; Wei et al., 2018; Zhou et al., 2018), authors determined δ_ET values by analogy to Ohm's law. The so-called “flux gradient method” (Lee et al., 2007) is based on the premise that the ET flux density rate (F_ET, ${\mathrm{L}}^{\mathrm{3}}{\mathrm{L}}^{-\mathrm{2}}{\mathrm{T}}^{-\mathrm{1}}$, expressed typically in $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) is proportional to $\mathrm{\Delta }{\mathit{\chi }}_{\text{atm}}/\mathrm{\Delta }{z}_{\text{atm}}$ (L⁻¹, typically in m⁻¹), the gradient of water vapour mixing ratio between two observation heights (with z_atm standing for height):

$$\begin{array}{}\text{(8)}& F_{ET}=-K{\displaystyle \frac{{\mathit{\rho }}_{\text{atm}}}{M_{\text{atm}}}}{\displaystyle \frac{\mathrm{\Delta }{\mathit{\chi }}_{\text{atm}}}{\mathrm{\Delta }{z}_{\text{atm}}}}.\end{array}$$

The water vapour transport is determined by the overall conductance of the air boundary layer expressed here as $-K{\mathit{\rho }}_{\text{atm}}/M_{\text{atm}}$ with ρ_atm (M L⁻³) and M_atm (M L⁻³, units of kg mol⁻¹) being the volumetric mass and molecular weight of dry air and K (L² T⁻¹) being the eddy diffusivity of water vapour. The isotopic ratio of ET (R_ET, –), which can be defined as the ratio of the flux density rates of the rare (superscript i) and abundant (superscript j) isotopologues (ⁱF_ET and ^jF_ET, respectively), can be therefore expressed as

$$\begin{array}{}\text{(9)}& R_{ET}={}^{\mathrm{i}}{F}_{ET}/{}^{\mathrm{j}}{F}_{ET}\approx \mathrm{\Delta }{}^{\mathrm{i}}{\mathit{\chi }}_{\text{atm}}/\mathrm{\Delta }{}^{\mathrm{j}}{\mathit{\chi }}_{\text{atm}},\end{array}$$

assuming that differences in K among water stable isotopologues are not significant, i.e. ${}^{\mathrm{i}}K={}^{\mathrm{j}}K=K$ (Yakir and Wang, 1996; Griffis et al., 2005). ⁱχ_atm and ^jχ_atm are the water vapour mixing ratio of rare and abundant isotopologues, respectively. Equation (9) can be further rearranged as

$$\begin{array}{}\text{(10)}& {}^{\mathrm{i}}{\mathit{\chi }}_{\text{atm}}=R_{ET}{}^{\mathrm{j}}{\mathit{\chi }}_{\text{atm}}+C,\end{array}$$

where R_ET is the slope and C (–) is the y intercept of the linear relationship between ⁱχ_atm and ^jχ_atm. Equation (10) becomes in δ notation

$$\begin{array}{}\text{(11)}& {\mathit{\delta }}_{\text{atm}}={\mathit{\delta }}_{ET}+C/R_{\text{std}}{\displaystyle \frac{\mathrm{1}}{{}^{\mathrm{j}}{\mathit{\chi }}_{\text{atm}}}}\end{array}$$

by dividing its left and right terms by ^jχ_atmR_std with R_std, the isotopic ratio of the internationally accepted water standard, namely the Vienna Standard Mean Ocean Water (V-SMOW) (Gonfiantini, 1978). We note that, by assuming ^jχ_atm≈χ_atm, the flux gradient and Keeling plot techniques are mathematically identical if $C={\mathit{\chi }}_{\text{bg}}({\mathit{\delta }}_{\text{bg}}-{\mathit{\delta }}_{ET})R_{\text{std}}=s\cdot R_{\text{std}}$, with s being the Keeling plot slope.

Griffis et al. (2010) and Good et al. (2012) used a combination of the EC technique and infrared tunable-diode-laser (TDL) water isotope spectroscopy to derive δ_ET values from simultaneous changes in wind velocity (ω, L T⁻¹) and ⁱχ_atm. In this statistical framework and by

• i.

  considering that air density and storage fluctuations are negligible during the measurement period (typically 30 min) and

• ii.

  changing the coordinate system so that the vertical wind velocity mean value ($\stackrel{\mathrm{‾}}{\mathit{\omega }}$) equals zero,

F_ET is expressed as

$$\begin{array}{}\text{(12)}& F_{ET}={\displaystyle \frac{{\mathit{\rho }}_{\text{atm}}}{M_{\text{atm}}}}\stackrel{\mathrm{‾}}{{\mathit{\omega }}^{\prime }{\mathit{\chi }}_{\text{atm}}^{\prime }}.\end{array}$$

The term $\stackrel{\mathrm{‾}}{{\mathit{\omega }}^{\prime }{\mathit{\chi }}_{\text{atm}}^{\prime }}$ is the average (overbar symbol) product of the differences between instantaneous and mean values (indicated by the prime symbols) of wind velocity and water vapour mixing ratio, in other words the covariance between the ω and χ_atm monitored variables. Similar to Eq. (11), we obtain after converting into δ notation the expression for the isotopic composition of ET:

$$\begin{array}{}\text{(13)}& {\mathit{\delta }}_{ET}={\displaystyle \frac{{}^{\mathrm{i}}{F}_{ET}/{}^{\mathrm{j}}{F}_{ET}}{R_{\text{std}}}}-\mathrm{1}={\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{\omega }}^{\prime }{}^{\mathrm{i}}{\mathit{\chi }}_{\text{atm}}^{\prime }}/\stackrel{\mathrm{‾}}{{\mathit{\omega }}^{\prime }{}^{\mathrm{j}}{\mathit{\chi }}_{\text{atm}}^{\prime }}}{R_{\text{std}}}}-\mathrm{1}.\end{array}$$

An alternative to Eq. (13) consists of considering the high-frequency variations of δ_atm rather than those of the individual mixing ratios ⁱχ_atm and ^jχ_atm. For this the isoflux (Lee et al., 2009), defined as $\stackrel{\mathrm{‾}}{{\mathit{\omega }}^{\prime }{\mathit{\delta }}_{\text{atm}}^{\prime }}$ (${\mathrm{L}}^{\mathrm{3}}{\mathrm{L}}^{-\mathrm{2}}{\mathrm{T}}^{-\mathrm{1}}$), is introduced:

$$\begin{array}{}\text{(14)}& {\mathit{\delta }}_{ET}={\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{\text{atm}}}}{\stackrel{\mathrm{‾}}{{\mathit{\omega }}^{\prime }{\mathit{\chi }}_{\text{atm}}^{\prime }}}}\stackrel{\mathrm{‾}}{{\mathit{\omega }}^{\prime }{\mathit{\delta }}_{\text{atm}}^{\prime }}+\stackrel{\mathrm{‾}}{{\mathit{\delta }}_{\text{atm}}}.\end{array}$$

3.1.2 Progress and challenges

In a review of isotope techniques for determination of the concomitant flux and isotopic composition of evapotranspiration, Griffis (2013) summarized the inherent limitations of the Keeling plot technique from the literature. The general assumption that atmospheric water vapour and its isotopic composition result from the turbulent mixing of only two sources was reported to be often violated. Reasons for this may be strong vertical gradients of the water vapour mixing ratio and isotopic composition or strong differences between δ_E and δ_E leading to the emergence of diffusion and air entrainment processes (Lee et al., 2006, 2012). The spatial Keeling plot suffers particularly from the fact that the different heights at which δ_atm is measured are representative for differently footprints areas of the studied ecosystem. While this may not be a problem for a homogeneous cropland, the reliability of the Keeling plot should be generally questioned for mixed vegetation (such as represented in Fig. 1) with strong lateral variabilities not only in δ_atm and χ_atm but also in soil water isotopic composition. In addition, the application of the spatial Keeling plot should not be conditioned based on a wide span of χ_atm values only but naturally on the quality of its linear fit. Griffis (2013) argued as well that the flux gradient approach suffers from a narrow range of application; e.g. it may not be suitable in certain cases, such as below forest canopies or above tall vegetation.

Regardless of these limitations or complications, Good et al. (2012) and Hu et al. (2020) provided comprehensive comparisons of the various techniques (Keeling plot, flux gradient, and EC) for determination of δ_ET. In addition to the temporal and spatial variations of the Keeling plot, Good et al. (2012) tested a third option where, instead of instantaneous measurements of δ_atm and χ_atm collected during 30 min, the mean values of δ_atm and χ_atm are calculated at each observation height (n=4) and used for regression. After a detailed uncertainty analysis, they concluded that the use of mean values instead of individual data points increased the uncertainty associated with δ_ET, regardless of the kind (temporal vs. spatial) of Keeling plot. In addition, the techniques of the temporal and spatial Keeling plot yielded significantly different values of δ_ET for the same time interval. The authors could not conclude which value was the most representative. In addition, they found a good agreement between the Keeling plot technique, applied at different heights, and the flux gradient method due to the aforementioned mathematical similarities. Hu et al. (2020) compared at one irrigated maize crop δ_ET values determined with either the Keeling plot or the flux gradient approaches. They tested different regression models with the Keeling plot method, i.e. ordinary least-squares regression, geometric-mean regression, and York's solution (for details see Pataki et al., 2003; Wehr and Saleska, 2017). These models differ in the way errors made on $\mathrm{1}/{\mathit{\chi }}_{\text{atm}}$ and δ_atm (see Eq. 5) relate to each other and whether they may be considered constant over their measurement ranges. As such, they yield differences in δ_ET estimates. Hu et al. (2020) could illustrate the necessity of choosing an appropriate regression model that reflects the dependency of spectrometer-specific errors on the water vapour mixing ratio. Yepez et al. (2005) and Wang et al. (2013) combined the Keeling plot technique with their closed chamber systems. During the course of measurement (e.g. 6 min in the study of Yepez et al., 2005) and for the Keeling plot approach to be valid, the increase of the chamber water vapour mixing ratio (10–15 mmol mol⁻¹ in Yepez et al., 2005) should not induce changes in both ET flow rate and isotopic composition. The fulfilment of this requirement of the Keeling plot technique is verified in a first approach by the very existence of a linear relationship between chamber air $\mathrm{1}/\mathit{\chi }$ and δ values. However, it could be argued that the linear form of the regression equation should survive a linear change in either $\mathrm{1}/{\mathit{\chi }}_{\text{atm}}$ or δ_atm. Another issue related to the use of chamber systems is the occurrence of water vapour condensation on the inside of the chamber or within the tubing system, e.g. following changes of incoming solar radiation during measurement. This may result in eventual isotopic fractionation leading to unreliable (i.e. unstable and underestimated) observations of chamber air δ values. To avoid such problems, the volume of the chamber is critical (i.e. the bigger it is, the less sensitive it is to abrupt changes of outside conditions), and active ventilation is mandatory. Ventilation not only prevents from condensation problems and pressure anomalies (Longdoz et al., 2000) but also guarantees the prevalence of turbulent mixing conditions in the chamber. The latter may not be ensured by a high turnover rate alone, i.e. the ratio of chamber volume and flow rate of flushed air. It is an important prerequisite of the application of both techniques based on the Keeling plot and alternating measurements of the water vapour mixing ratio and isotopic composition of inlet and outlet air (Eq. 7). Measurements with dynamically purged chambers, which are combined with the latter type of mass balance applications, may reduce the problem of condensation inside the chamber. A possibility is to flush the chamber with dry air so that the increase in water vapour mixing ratio and (positive or negative) change of the isotopic composition measured at the outlet relative to the inlet directly reflect the volume and isotopic composition of the moisture added by ET. Stable measurements over a certain time period, depending on both chamber volume and inflow rate, would indicate ISS, and δ_ET may be directly measured without any further calculations (e.g. Wang et al., 2010). However, dry air can stress the enclosed plants by artificially increasing the chamber air vapour pressure deficit, which ultimately can result in NSS conditions. In this case, a steady increase of chamber air χ should not be observed during the course of the measurement, as it would be a sign of a significant difference of micrometeorological conditions (temperature, vapour pressure deficit, and wind speed values) between ambient and chamber air.

As stated above, both the techniques of the temporal Keeling plot and the flux gradient suffer from the need of a high spatial gradient in the water vapour mixing ratio and isotopic composition between the soil or canopy surface and the free atmosphere to obtain precise values of δ_ET. One alternative to sampling water vapour in atmospheric profiles at fixed heights is to use a small (about a few metres high) field lift system, the modus operandi of which is based on the principles established by Mayer et al. (2009) and Noone et al. (2013), for a continuous monitoring of atmospheric height profiles of the water vapour isotopic composition. To the authors' knowledge, only one study on an evergreen forest made use of the principle in the context of ET partitioning (Berkelhammer et al., 2016). Ney and Graf (2018) designed a portable lift system for measuring the atmospheric water vapour and CO₂ mixing ratios in the field for various crops at a half-hourly temporal resolution. Their system should allow for measuring highly vertically resolved water vapour isotopic profiles. For this, however, high-throughput and high-frequency isotopic analysers are needed to provide reliable information on ecosystem fluxes. Commercially available cavity ring-down laser spectrometers operate at low flow rate (φ) and frequency (f) (e.g. $\mathrm{25}\le \mathit{\phi }\le$ 35 mL min⁻¹ and f≈ 1 Hz for the L2120-i, L2130-i, and L2140-i by Picarro, Inc., Santa Clara, CA, USA) and are, thus, not suitable for such measurements. Other isotope analysers, such as the quantum-cascade-laser (QCL) trace gas monitor (Aerodyne, Inc., Billerica, MA, USA; φ≤ 250 L min⁻¹ and f = 10 Hz) have the potential to fulfil the abovementioned requirements.

Compared to the Keeling plot and flux gradient approaches, the eddy covariance technique derives from micrometeorological theory (first principles). Where applicable, this makes it a solid alternative less subjected to assumptions. However, as a result of its high data acquisition rate and associated noise, the EC technique provides δ_ET estimates with higher uncertainty, largely determined by random measurement errors (Hollinger and Richardson, 2005; Loescher et al., 2006; Rannik et al., 2016). Good et al. (2012) determined this uncertainty to be proportional to the inverse of the correlation coefficient between ω and χ_atm, i.e. the covariance of ω and χ_atm divided by the product of their measurement errors.

One important feature of the EC isotope technique resides in its ability to provide δ_ET estimates at the field scale and therefore demarks itself from the abovementioned approaches. Griffis et al. (2010) and Griffis et al. (2011) demonstrated the reliability of δ_ET data obtained with the eddy covariance technique from the agreement between measurements of the water vapour mixing ratio and ET flux of their traditional infrared analyser (LI-7500, LI-COR, Inc., Lincoln, NE, USA) and a fast-response and high-frequency water isotope spectrometer, i.e. the lead–salt tunable-diode-laser spectrometer TGA200A (Campbell Scientific, Inc., Logan, UT, USA; φ = 1.7 L min⁻¹ at f = 10 Hz). However, to our knowledge, the production of this instrument has been discontinued. Recently, Braden-Behrens et al. (2019, 2020) showed that EC measurements could be performed using a high-flow (φ≈ 4.2 L min⁻¹) laser spectrometer clocked at 2 Hz only (2 Hz-HF-WVIA, Los Gatos Research Inc., San Jose, CA, USA). They underlined the importance of heating the sampling tubing at the point of intake in order to avoid problems of condensation and high-frequency dampening as shown by spectra and cospectra.

3.2 Isotopic composition of evaporation

3.2.1 Methods

Two options are found in the literature (Fig. 2i) for determining the isotopic composition of the E flux, δ_E:

• i.

  by solving one of either mass balance equations (Eqs. 7 or 11; see Sect. 3.1) in combination with dynamically purged closed bare-soil chambers (15 % of the reviewed studies, e.g. Dubbert et al., 2013, 2014b) or

• ii.

  by solving the so-called “Craig and Gordon equation” (Eq. 18 below), which is derived from the atmospheric part of a transport model of water stable isotopologues, based on an analogy to Ohm's law (Craig and Gordon, 1965) (69 % of the studies).

The two approaches differ in numerous aspects: while the first is non-destructive and requires online and continuous measurements of a few variables (i.e. water vapour mixing ratio and isotopic composition of the chamber inlet and outlet air), the second relies – with the exception of the study of Quade et al. (2019) – on destructive sampling of the soil and offline analysis of the extracted water. The Craig and Gordon equation demarks itself from Eqs. (7) and (11) also due to its complex parametrization. Craig and Gordon (1965) classically interpreted the temporal changes in δ_E of a free water body with the help of a linear-resistance model. We will shortly present the widely used model variation for water bound to the soil media (for an in-depth review, the reader is kindly referred to Horita et al., 2008). The only significant difference to the original model is the evaporating front vertical coordinate (z_EF), which may not correspond to that of the soil surface depending on the evaporation stage (Or et al., 2013; Merz et al., 2018). The isotopic ratio of evaporation, R_E, is expressed as the ratio of ⁱF_E and ^jF_E, i.e. the water vapour flux density rates (${\mathrm{L}}^{\mathrm{3}}{\mathrm{L}}^{-\mathrm{2}}{\mathrm{T}}^{-\mathrm{1}}$) in rare and abundant isotopologues, respectively, originating from the EF:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{R}_{\mathrm{E}}={\displaystyle \frac{{}^{\mathrm{j}}{F}_{\mathrm{E}}}{{}^{\mathrm{j}}{F}_{\mathrm{E}}}}& ={\displaystyle \frac{\mathrm{\Delta }{}^{\mathrm{i}}{\mathit{\chi }}_{\text{atm}}}{\mathrm{\Delta }{}^{\mathrm{j}}{\mathit{\chi }}_{\text{atm}}}}\\ & ={\displaystyle \frac{\mathrm{1}}{{}^{\mathrm{i}}r/{}^{\mathrm{j}}r}}\cdot {\displaystyle \frac{{}^{\mathrm{i}}{\mathit{\chi }}_{\text{atm}}\left(z_{\text{EF}}\right)-{}^{\mathrm{i}}{\mathit{\chi }}_{\text{atm}}\left(z_{\text{atm}}\right)}{{}^{\mathrm{j}}{\mathit{\chi }}_{\text{atm}}\left(z_{\text{EF}}\right)-{}^{\mathrm{j}}{\mathit{\chi }}_{\text{atm}}\left(z_{\text{atm}}\right)}}.\end{array}\end{array}$$

We note that Eq. (15) is analogous to Eq. (9) (Lee et al., 2007; see Sect. 3.1), with the exception that the bulk resistances to vapour transport of the rare and abundant isotopologues (ⁱr and ^jr, respectively, T L⁻¹) are not assumed equal. It follows from the fact that ⁱr and ^jr relate to the air layer delimited between z_EF and z_atm (and not between the two observation heights in Eq. 9), where not only purely turbulent transport or eddy diffusivity but also molecular diffusion and laminar flow are relevant. Furthermore, Craig and Gordon (1965) conceptualized the existence of a water-vapour-saturated (superscript “sat”) air layer at the EF where isotopic thermodynamic equilibrium prevails:

$$\begin{array}{}\text{(16a)}& {}^{\mathrm{j}}{\mathit{\chi }}_{\text{atm}}\left(z_{\text{EF}}\right)={}^{\mathrm{j}}{{\mathit{\chi }}_{\text{atm}}}^{\text{sat}},\text{(16b)}& {}^{\mathrm{i}}{\mathit{\chi }}_{\text{atm}}\left(z_{\text{EF}}\right)={}^{\mathrm{i}}{{\mathit{\chi }}_{\text{atm}}}^{\text{sat}}={}^{\mathrm{j}}{{\mathit{\chi }}_{\text{atm}}}^{\text{sat}}{R}^{\text{sat}}={}^{\mathrm{j}}{{\mathit{\chi }}_{\text{atm}}}^{\text{sat}}{\displaystyle \frac{R_{\text{EF}}}{{\mathit{\alpha }}_{\text{eq}}}}.\end{array}$$

where R^sat and R_EF are the isotopic ratios of the saturated air layer and of the soil liquid water at the EF, respectively. α_eq (–) is the isotopic equilibrium fractionation factor, first empirically determined by Horita and Wesolowski (1994). α_eq depends on the soil temperature at the EF (T_EF, K):

$$\begin{array}{}\text{(17)}& {\mathit{\alpha }}_{\text{eq}}\left(T_{\text{EF}}\right)=\exp \left({\displaystyle \frac{A}{T_{\text{EF}}^{\mathrm{2}}}}+{\displaystyle \frac{B}{T_{\text{EF}}}}+C\right),\end{array}$$

with constants $A=\mathrm{24.844}\times {\mathrm{10}}^{\mathrm{6}}$, $B=-\mathrm{76.248}\times {\mathrm{10}}^{\mathrm{3}}$, and C = 0.052612 for ²H and $A=\mathrm{1.137}\times {\mathrm{10}}^{\mathrm{6}}$, $B=-\mathrm{0.4156}\times {\mathrm{10}}^{\mathrm{3}}$, and C = −0.0020667 for ¹⁸O. Craig and Gordon (1965) identified the ratio of bulk resistances ${}^{\mathrm{i}}r/{}^{\mathrm{j}}r$ as the isotopic kinetic fractionation factor (α_K, –). Finally, by

• i.

  considering that ⁱχ_atm(z_atm)=^jχ_atm(z_atm)R_atm,

• ii.

  dividing the right-hand term numerator and denominator of Eq. (15) by ${}^{\mathrm{j}}{{\mathit{\chi }}_{\text{atm}}}^{\text{sat}}$, and

• iii.

  converting R_E into δ_E, we obtain:

$$\begin{array}{}\text{(18)}& {\mathit{\delta }}_E={\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{K}}(\mathrm{1}-h_{\text{atm}})}}\left({\displaystyle \frac{{\mathit{\delta }}_{\text{EF}}+\mathrm{1}}{{\mathit{\alpha }}_{\text{eq}}}}-({\mathit{\delta }}_{\text{atm}}+\mathrm{1})h\right)-\mathrm{1},\end{array}$$

where h_atm (–) is the relative humidity of the ambient atmosphere measured at vertical coordinate z_atm and defined as the ratio ${}^{\mathrm{j}}{\mathit{\chi }}_{\text{atm}}\left(z_{\text{atm}}\right)/{}^{\mathrm{j}}{{\mathit{\chi }}_{\text{atm}}}^{\text{sat}}$. The possible difference in temperature measured at z_atm and z_EF should be accounted for by normalizing h_atm to the saturated vapour pressure ($\mathrm{M}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{T}}^{-\mathrm{2}}$, usually expressed in Pa) at the temperature of the EF (Rothfuss et al., 2010; Quade et al., 2019).

Craig and Gordon (1965) argued that the kinetic fractionation factor was inversely proportional to the ratio of the molecular diffusivities of ${}^{\mathrm{1}}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{16}}\mathrm{O}$ (^jD) and of either ¹H²H¹⁶O or ${}^{\mathrm{1}}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$ (ⁱD):

$$\begin{array}{}\text{(19)}& {\mathit{\alpha }}_{\mathrm{K}}={\displaystyle \frac{{}^{\mathrm{i}}r}{{}^{\mathrm{j}}r}}={\displaystyle \frac{{}^{\mathrm{j}}D}{{}^{\mathrm{i}}D}}.\end{array}$$

Figure 5 Effect of the water status of the soil, i.e. the positioning of the evaporating front (EF, dashed line), on the value of the kinetic fractionation factor (α_K). Panels (a–c) refer to the situation of a saturated soil (subscript “wet”) where the EF is located at the soil surface; panels (d–f) refer to a dry soil with the EF below the soil surface. The corresponding soil water total (tot; solid line), liquid (liq; dotted line), and vapour (vap; dash-dotted line) isotopic flux profiles (E_i, M L⁻³) (a, f); soil liquid isotopic composition profiles (b, e) are reported as well. Adapted from Braud et al. (2005).

Merlivat (1978) and later Luz et al. (2009) quantified the diffusivity ratios at 1.0251 and 1.0285 for ¹H²H¹⁶O or ${}^{\mathrm{1}}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$ isotopologues, respectively. Dongmann et al. (1974) (but see also Brutsaert, 1975) extended Eq. (19) to different aerodynamic regimes in the air boundary layer delimited by z_EF and z_atm:

$$\begin{array}{}\text{(19')}& {\mathit{\alpha }}_{\mathrm{K}}={\left({\displaystyle \frac{{}^{\mathrm{j}}D}{{}^{\mathrm{i}}D}}\right)}^n,\end{array}$$

where n (–) is an unitless factor ranging from 0.5 (corresponding to fully turbulent conditions) to 1 (fully diffusive), with a value of $\mathrm{2}/\mathrm{3}$ representative of laminar-flow conditions. Mathieu and Bariac (1996) proposed to define n in the case of evaporation from soil as a linear function of soil volumetric water content observed in the surface layer (θ_surf; L³ L⁻³, typically in cm³ cm⁻³). n would range between 0.5 when θ_surf reaches θ_sat, the water content value at saturation, and 1 for θ_surf=θ_res, the value of residual water content (see Fig. 5).

$$\begin{array}{}\text{(20)}& n={\displaystyle \frac{({\mathit{\theta }}_{\text{surf}}-{\mathit{\theta }}_{\text{res}})n_{\text{atm}}+({\mathit{\theta }}_{\text{sat}}-{\mathit{\theta }}_{\text{surf}})n_{\text{soil}}}{{\mathit{\theta }}_{\text{sat}}-{\mathit{\theta }}_{\text{res}}}}\end{array}$$

In Mathieu and Bariac's conceptual framework the establishment of a dry soil surface layer results in added isotopic resistance by increasing the relative importance of gaseous molecular diffusion (i.e. in the tortuous soil pore network) in the overall transport of water vapour from the EF towards the well-mixed atmosphere. In the case of a fully water-saturated layer in contact with the free atmosphere, the opposite happens: water vapour leaving the rough surface is preferentially transported in a turbulent manner, leading to smaller n values.

3.2.2 Progress and challenges

To calculate δ_E with the Craig and Gordon equation requires simultaneous observations of h_atm, T_EF, δ_atm, and δ_EF. The first two variables are typically monitored with e.g. capacitive sensing. As for δ_atm, its value is determined from online or offline isotopic analysis after sampling of the atmospheric water vapour (see Sect. 3.1).

The variable most challenging to estimate is δ_EF (Fig. 5b and e). It greatly depends on how soil is sampled at the EF. However, there is no consensus on how this should be done in the literature (see the column on isotopic measurements in Table A1). Some studies do not precisely report the soil depth, which is considered to be the EF (e.g. Wang and Yakir, 2000; Yepez et al., 2003; Williams et al., 2004). In other studies (Yepez et al., 2005; Zhang et al., 2011; Dubbert et al., 2013) soil profiles are partially or entirely sampled at higher vertical (centimetre) resolution. Pioneer works on isotopic transport in saturated or non-saturated isothermal soils under steady-state evaporation (Zimmermann et al., 1967; Allison, 1982; Barnes and Allison, 1983) showed that the EF, i.e. the theoretical and continuous boundary between the soil “regions” dominated by either liquid or vapour flow (Fig. 5a and f), is associated with the highest isotopic composition (${\mathit{\delta }}_{\text{soil}}^{\text{liq}}$) value of the liquid soil water (Fig. 5d–f). Later this family of models was extended to unsaturated soil water conditions, non-isothermal conditions, and time-variable evaporation flux (e.g. Barnes and Allison, 1988; Barnes and Walker, 1989). More recently, Braud et al. (2005) and Haverd and Cuntz (2010) implemented isotopic transport in both liquid and vapour phases of the soil, with a coupling to temperature dynamics, in numerically solved SVAT models (SiSPAT-Isotope and Soil–Litter–Iso). All the abovementioned studies underline the localized character of the EF and the strong isotopic gradient in liquid water at its location. The determination of the EF location may be problematic, especially in the case of a receding EF (z_EF≠z_surf, Fig. 5d), which is generally the case in arid regions between rare precipitation events. Thus, sampling soil roughly from the surface does not allow for a precise determination of δ_EF and may lead to errors in δ_E estimates. Rothfuss et al. (2010) could demonstrate for a well-watered soil (i.e. with z_EF=z_surf, Fig. 5b) that sampling of only a few centimetres of soil at the surface and using the corresponding δ_surf in Eq. (18) could lead to a significant underestimation of δ_E. This would lead in turn to an overestimation of $T/ET$, since negative changes in δ_E translate into positive changes in $T/ET$, i.e.$\frac{\partial \left(\frac{T}{ET}\right)}{\partial \left({\mathit{\delta }}_E\right)}=\frac{{\mathit{\delta }}_{ET}-{\mathit{\delta }}_{\mathrm{T}}}{({\mathit{\delta }}_{\mathrm{T}}-{\mathit{\delta }}_E{)}^{\mathrm{2}}}<\mathrm{0}$ (when δ_ET < δ_E, which is generally the case). The spatial (vertical) resolution of the soil sampling should therefore be as high as possible to be able to identify z_EF precisely. For their specific case, Brunel et al. (1997) estimated also that the determination of the δ_EF value was the greatest source of uncertainty of $T/ET$.

After sampling in the field, water is recovered from the soil in the laboratory using one of six extraction methods: cryogenic vacuum extraction (Araguás-Araguás et al., 1995; West et al., 2006), azeotropic distillation (Revesz and Woods, 1990), direct vapour equilibration (Wassenaar et al., 2008), high-pressure mechanical squeezing (Kelln et al., 2001), centrifugation (Mubarak and Olsen, 1976), or microwave extraction (Munksgaard et al., 2014). Other methods include the use of soil liquid water samplers (Wenninger et al., 2010; Sutanto et al., 2012). Finally, δ_EF is measured by isotope ratio mass spectrometry (IRMS) or isotope ratio infrared spectrometry (IRIS). Note that an alternative consists of letting soil water directly equilibrate with CO₂ without the need for water extraction (one study is Ferretti et al., 2003, which follows the method of Scrimgeour, 1995). In this framework, pure CO₂ is injected in the exetainer containing the soil sample following evacuation. After a 3 d long water–CO₂ equilibration period, the δ¹⁸O value of CO₂ is measured by isotope mass spectrometry and used to infer that of water at equilibrium. Orlowski et al. (2016a, b) provided evidence from laboratory benchmarks of the different techniques that the isotopic composition of the recovered water could be sensitive to the extraction approach and extraction time as well as to the soil type and values of water and organic content. The same authors also observed that IRMS and IRIS techniques yielded different results in general and especially for clay loam soil water, which they related to interferences in the absorption spectrum during analysis with the latter technique. In addition, Orlowski et al. (2018) concluded from a worldwide inter-comparison of cryogenic-vacuum-extraction facilities that the general consensus in the isotopic ecohydrology community, stating that cryogenic vacuum extraction is the standard water recovery technique, should be questioned. Orlowski et al. (2016a, b, 2018) highlighted the limitations of the most popular extraction approach, i.e. based on the combination of destructive sampling and vacuum extraction (69 % of the reviewed studies), which calls for the development of other techniques for a precise quantification of δ_EF.

In the last few years Rothfuss et al. (2013), Volkmann and Weiler (2014), and Gaj et al. (2016) successfully validated and tested alternatives to destructive sampling and offline isotopic analysis approaches. They developed systems based on the combination of gas-permeable membranes (e.g. rigid hydrophobic microporous polypropylene, Membrana GmbH, Germany, or polyethylene, Porex Technologies, Aachen, Germany) with laser-based spectrometry for the non-destructive collection of the soil atmosphere and the online monitoring of its water vapour isotopic composition (${\mathit{\delta }}_{\text{soil}}^{\text{vap}}$). For this, the soil atmosphere is either

• i.

  flushed with a carrier gas (dry synthetic air, i.e. 20.5 % in N₂, or 100 % N₂) at low flow rate in the range of 50–100 mL min⁻¹ (Rothfuss et al., 2013; Volkmann and Weiler, 2014; Gaj et al., 2016) or

• ii.

  extracted with a vacuum pump (Volkmann and Weiler, 2014).

Both modi operandi allow for long-term and repeated measurements across the soil profile provided that condensation is avoided in the sampling line. For this, the collected air, which is (quasi-)saturated with water vapour, is diluted with the carrier gas and the sampling lines are heated, if necessary (Quade et al., 2019; Kühnhammer et al., 2019). Rothfuss et al. (2013) observed near-isotopic-equilibrium conditions between liquid and vapour in the soil pore space and provided temperature calibration equations yielding results analogous to those of Horita and Wesolowski (1994) for converting ${\mathit{\delta }}_{\text{soil}}^{\text{vap}}$ into ${\mathit{\delta }}_{\text{soil}}^{\text{liq}}$ values. They further show that isotopic equilibrium conditions still prevailed at low soil volumetric water content, possibly also for soil water vapour relative humidity values lower than 1. Their method was successfully applied to laboratory experiments with sand (Gangi et al., 2015; Rothfuss et al., 2015) and silt loam (Quade et al., 2018). Oerter et al. (2017) compared ${\mathit{\delta }}_{\text{soil}}^{\text{liq}}$ values estimated with the monitoring method of Rothfuss et al. (2013) on the one hand and the direct equilibrium and vacuum extraction methods on the other hand. They found a good correlation between the two approaches (root mean square error – RMSE – equal to 0.6 ‰ for δ¹⁸O and within 1.7 ‰–3.1 ‰ for δ²H). Volkmann and Weiler (2014) tested their own design of a water vapour probe under field conditions and could show that it produced ${\mathit{\delta }}_{\text{soil}}^{\text{liq}}$ values in agreement with those following destructive sampling and isotopic analysis with the direct equilibration method (Garvelmann et al., 2012). The inter-method (destructive vs. non-destructive) RMSE values were comparable to the intra-method variability of soil water δ values. The latter variability could not be disentangled into systematic methodological error and natural (lateral) heterogeneity in soil water isotopic composition. Kübert et al. (2020) conducted a comparison study of the method of Rothfuss et al. (2013) with cryogenic vacuum extraction and centrifugation during an irrigation pulse-labelling experiment in a semi-natural temperate grassland. They highlighted that the non-destructive method could capture temporal dynamics of the isotopic composition, while destructive sampling included both the temporal change and spatial heterogeneity.

To date there are two ET partitioning studies, in which δ_E was determined from non-destructive isotopic analysis using soil-liquid-water–water-vapour equilibration. Quade et al. (2019) applied the method of Rothfuss et al. (2013) in a sugar beet field in a temperate climate (Germany), while Gaj et al. (2016) used commercially available soil gas probes (BGL-30, METER Group, Munich, Germany), following the same modus operandi as Volkmann and Weiler (2014), during a field study in central Namibia under semi-arid conditions. Such applications are promising for the specific purpose of partitioning ET, as they provide insights into subdaily dynamics of δ_E from the online assessment of the positioning and isotopic composition of water at the EF. However, one noticeable disadvantage is the need for deploying a laser spectrometer at the experimental site. A possible way around has been recently proposed by Havranek et al. (2020) as a compromise: water vapour samples are collected and stored automatically in flasks from the soil profile in the field following the approach of Rothfuss et al. (2013) and transported back to the laboratory where the isotopic analyses are performed.

Table 1 Effect of the consideration of non-saturated soil water vapour phase on the estimation of the isotopic composition of evaporation (δ_E) using the model of Craig and Gordon (1965). Conditions of the transport of pure diffusive water vapour (n=1) prevail, leading to values of the kinetic fractionation factor (α_K) of 1.0251 und 1.0285 for ²H and ¹⁸O. Values for T_EF, h_atm, δ_EF, and δ_atm are chosen exemplarily.

Another important factor that influences the precision of δ_E estimates is the choice of the value of the kinetic fractionation factor, α_K. Only a handful of studies attempted to estimate or model α_K for soil E. Braud et al. (2009) simulated α_K values during long-term laboratory experiments with the SVAT model SiSPAT-Isotope. They found a decreasing trend in α_K value from saturated to unsaturated soil conditions, which contradicts the model of Mathieu and Bariac (1996). Results similar to the study by Braud et al. (2009) were obtained by Rothfuss et al. (2015) during a long-term soil column laboratory experiment. Quade et al. (2018) tested two different methods for quantifying α_K during a series of bare-soil evaporation experiments on monoliths (100 L soil volume) under semi-controlled conditions, i.e. the following:

• i.

  The first method is the inversion of the Craig and Gordon equation (Eq. 18) in a single isotope-framework (i.e. based on either δ¹⁸O or δ²H values) with input variables h_atm, T_EF, δ_atm, δ_EF, and δ_E.

• ii.

  The second method is the inversion of the Craig and Gordon equation in a dual-isotope framework. More specifically, α_K is determined from the approximation of the slope (S, –) of the soil E line (S_E = $\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{2}}\mathrm{H}}_{\text{soil}}^{\text{liq}}/\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{soil}}^{\text{liq}}$, –) in a [δ¹⁸O, δ²H] coordinate system following Gat (2000):

  $$\begin{array}{}\text{(21)}& S_{\mathrm{E}}\left(t\right)={\displaystyle \frac{{\left[h\left(t\right)\left({\mathit{\delta }}_{\text{atm}}\left(t\right)-{\mathit{\delta }}_{\text{EF}}(t-\mathrm{1})\right)+{\mathit{\epsilon }}_{\text{eq}}\left(t\right)+\mathrm{\Delta }\mathit{\epsilon }\left(t\right)\right]}_{{}^{\mathrm{2}}\mathrm{H}}}{{\left[h\left(t\right)\left({\mathit{\delta }}_{\text{atm}}\left(t\right)-{\mathit{\delta }}_{\text{EF}}(t-\mathrm{1})\right)+{\mathit{\epsilon }}_{\text{eq}}\left(t\right)+\mathrm{\Delta }\mathit{\epsilon }\left(t\right)\right]}_{{}^{\mathrm{18}}\mathrm{O}}}},\end{array}$$

  with t being the time stamp and Δε (–) being the so-called kinetic effect, defined as

  $$\begin{array}{}\text{(22)}& \mathrm{\Delta }\mathit{\epsilon }=(\mathrm{1}-h_{\text{atm}})({}^{\mathrm{j}}D/{}^{\mathrm{i}}D-\mathrm{1})n,\end{array}$$

with superscripts i and j standing for the least and most abundant isotopologues, respectively. Equation (21) is solved in an implicit manner; in other words, S_E values simulated for time stamp t depend on δ_EF observation made at time stamp (t−1). The n value is then extracted from Eq. (21) from the confrontation of measured and simulated S_E, and finally fed into Eq. (19') to retrieve α_K values. Quade et al. (2018) showed that α_K could not be considered a constant value solely depending on flow conditions as proposed by Dongmann et al. (1974) or determined from soil water content following Mathieu and Bariac (1996) (Eq. 20). The second approach yielded α_K values in agreement with the literature (e.g. Merlivat, 1979). Quade et al. (2018) concluded that turbulent transport still played a significant role during the evaporation process, also under non-saturated conditions. These studies show that further sensitivity analyses of α_K to environmental conditions are needed to provide realistic estimates of δ_E and ultimately of $T/ET$ values. To our knowledge, there is no ET partitioning study in the field where α_K was considered to dynamically change (other than via the model of Mathieu and Bariac, 1996) depending on the contribution of air turbulence to water vapour transport in the free and canopy atmosphere, e.g. from measurements of the wind profile within and above the canopy (Brutsaert, 1975).

Another source of uncertainty arises from a lack of precise knowledge of the state of water vapour saturation at the EF (Braud et al., 2005; Rothfuss et al., 2015). In the Craig and Gordon equation, the kinetic fractionation factor is weighed by the term (h_EF − h_atm), where h_EF is generally assumed to be equal to 1, representative of saturated conditions at the EF. However, this assumption may not stand for dry soils considering the relationship between soil water matric potential ψ_EF ($\mathrm{M}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{T}}^{-\mathrm{2}}$, typically expressed in hectopascals or centimetres of water height) and pore space relative humidity at the EF (h_EF), as given by Kelvin's law:

$$\begin{array}{}\text{(23)}& h_{\text{EF}}=\exp \left({\displaystyle \frac{{\mathit{\psi }}_{\text{EF}}{M}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{w}}{R}_{\text{gas}}{T}_{\text{EF}}}}\right).\end{array}$$

M_w and ρ_w (M L⁻³) are the molar and volumetric masses of water, respectively, and R_gas ($\mathrm{M}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{T}}^{-\mathrm{3}}$) is the universal gas constant. Table 1 presents three different degrees of saturation of the soil vapour phase under isothermal conditions (T_EF = 20 ^∘C) and their corresponding hydrogen and oxygen isotopic composition values of the E flux (δ²H_E and δ¹⁸O_E). A decrease in h_EF from 100 % to 99.9 %, corresponding to an increase in the absolute ψ_EF value from 0 to 1000 hPa (i.e. from saturation to pF = 3) leads, for example, to an increase of 1.5 ‰ in δ²H_E and δ¹⁸O_E. A decrease in h_EF from 100 % to 99.3 % (increase from 0 to 10 000 hPa, i.e. pF = 4) would translate into an increase of 13 ‰ in δ²H_E and δ¹⁸O_E. Both δ²H_E and δ¹⁸O_E are affected in the same way by the change in value of the factor $\frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{K}}(h_{\text{EF}}-h_{\text{atm}})}$ (see Eq. 18), i.e. approximatively 2.0 ‰ per 0.1 % relative humidity. This may have a noticeable effect on the computation of $T/ET$, especially for δ¹⁸O, for which the difference of δ_E − δ_E is usually smaller than for δ²H. Gas exchange chambers and other experimental setups with semi-controlled conditions (such as weighing lysimeters) provide means to test the validity and existence of the abovementioned hypotheses and complications (e.g. Dubbert et al., 2013; Groh et al., 2018).

3.3 Isotopic composition of transpiration

3.3.1 Methods

The determination of the isotopic composition of T, δ_E, in the reviewed literature is mainly dependent on the underlying hypothesis of the isotopic steady or non-steady state (NSS) of T. While 42 % of all reviewed studies assume isotopic steady state (ISS), in other words that δ_E is time-invariant, 58 % do not make such an assumption but assume a transient state, i.e. NSS (Fig. 2j). This has substantial implications for the materials and methods used for the determination of δ_E. In the ISS approach, δ_E is directly inferred from the isotopic value of the leaf water source (δ_xyl), i.e. the water in the xylem vessels supplying the leaf water reservoir. This assumption is based on the fact that at ISS the flux density rate of the least abundant (ⁱF_xyl) (respectively most abundant, ^jF_xyl) isotopologue entering the leaf equals the flux density rate of the least abundant (ⁱF_T) (most abundant, ^jF_T) isotopologue leaving it by transpiration:

$$\begin{array}{}\text{(24a)}& {}^{\mathrm{j}}{F}_{\text{xyl}}={}^{\mathrm{j}}{F}_{\mathrm{T}},\text{(24b)}& {}^{\mathrm{i}}{F}_{\text{xyl}}={}^{\mathrm{i}}{F}_{\mathrm{T}}\Rightarrow {}^{\mathrm{j}}{F}_{\text{xyl}}{\mathit{\delta }}_{\text{xyl}}={}^{\mathrm{i}}{F}_{\mathrm{T}}{\mathit{\delta }}_{\mathrm{T}}\Rightarrow {\mathit{\delta }}_{\text{xyl}}={\mathit{\delta }}_{\mathrm{T}}.\end{array}$$

Note that in this framework an instantaneous change in ^jF_T, if compensated by a corresponding change in ^jF_xyl, should maintain the relationship δ_xyl=δ_T (Eq. 24b). In reality, a change in ^jF_T, due to variations in environmental factors (e.g. vapour pressure deficit of the free atmosphere and incoming solar radiation) implies a change in root water uptake depth profile, which in turn affects δ_xyl in the case of a heterogeneous distribution of the soil water isotopic composition (Rothfuss and Javaux, 2017). A new ISS is eventually reached, depending on the leaf water turnover time, i.e. the ratio of leaf water volume and transpiration rate (Dongmann et al., 1974; see below). To access xylem water, authors destructively sample stems (e.g. Wei et al., 2018; Quade et al., 2019), branches (e.g. Williams et al., 2004), or root water (Bijoor et al., 2011) and recover their water by e.g. cryogenic vacuum extraction.

The NSS approach for determining δ_E relies either on direct non-destructive monitoring (i.e. leaf chamber-based measurements, e.g. Wang et al., 2010) or on destructive sampling of plant material and subsequent extraction of water (e.g. Dubbert et al., 2013). In the former case, the modus operandi is the same as when operating ET and E chambers coupled to mass balance equations (see Sects. 3.1 and 3.2, respectively), except that one single leaf or several leaves are enclosed in the chamber (with a volume ranging from 150 to 190 cm³ in the literature), rather than the entire plant. It is then generally assumed that the leaf-scale δ_E estimate is also representative for the whole plant (e.g. Good et al., 2014). In the case of destructive sampling, δ_E is modelled on the basis of environmental factors (leaf temperature and free-atmosphere relative humidity) and isotopic variables. Two cases can be distinguished.

• i.

  δ_E is determined from the value of the isotopic composition of the leaf bulk water, δ_L, with the Craig and Gordon equation adapted to plant T (Sun et al., 2014; Hu et al., 2014):

  $$\begin{array}{}\text{(18')}& {\mathit{\delta }}_{\mathrm{T}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{K}}(\mathrm{1}-h)}}\left({\displaystyle \frac{{\mathit{\delta }}_{\mathrm{L}}+\mathrm{1}}{{\mathit{\alpha }}_{\text{eq}}}}-({\mathit{\delta }}_{\text{atm}}+\mathrm{1})h\right)-\mathrm{1}.\end{array}$$
• ii.

  The isotopic composition of leaf water may not be available, but that of its source, δ_xyl, is. The δ_E value is calculated after the relationship of Dongmann et al. (1974), which describes the temporal course of δ_L at a constant transpiration rate value (i.e. at permanent flow for T). The authors expressed the rate of change in δ_L as a function of the instantaneous difference between δ_xyl and δ_E at time t, by considering the leaf bulk water (delimited by volume per unit leaf area V_L, L³ L⁻²) to be transpired into ambient air at permanent flow (i.e. at density rate ^jF_T=^jF_xyl, as in Eq. 24a):

  $$\begin{array}{}\text{(25)}& \mathrm{d}{\mathit{\delta }}_{\mathrm{L}}={\displaystyle \frac{{}^{\mathrm{j}}{F}_{\mathrm{T}}}{V_{\mathrm{L}}}}\left({\mathit{\delta }}_{\text{xyl}}\right(t)-{\mathit{\delta }}_{\mathrm{T}}(t\left)\right)\mathrm{d}t.\end{array}$$

By combining Eqs. (18') and (25) and considering that δ_xyl is time-invariant, we obtain a first-order differential equation for δ_L, which yields after integration to

$$\begin{array}{}\text{(26)}& \begin{array}{rl}{\mathit{\delta }}_{\mathrm{L}}\left(t\right)=& {\mathit{\delta }}_{\mathrm{L}}(t\to +\mathrm{\infty })-({\mathit{\delta }}_{\text{xyl}}-{\mathit{\delta }}_{\mathrm{L}}(t\to +\mathrm{\infty }\left)\right)\\ & \cdot \exp \left(-{\displaystyle \frac{t}{{\mathit{\tau }}_{\mathrm{L}}}}{\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_{\text{eq}}{\mathit{\alpha }}_{\mathrm{K}}(\mathrm{1}-h_{\text{atm}})}}\right),\end{array}\end{array}$$

where the leaf water turnover time, τ_L, is defined as the ratio $\frac{V_{\mathrm{L}}}{{}^{\mathrm{j}}{F}_{\mathrm{T}}}$ and ${\mathit{\delta }}_{\mathrm{L}}(t\to +\mathrm{\infty })={\mathit{\delta }}_{\text{L\_ISS}}$, the isotopic composition of leaf bulk water when an isotopic steady state is reached. The latter term is expressed as

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{\mathit{\delta }}_{\text{L\_ISS}}=& {\mathit{\alpha }}_{\text{eq}}\left[{\mathit{\alpha }}_{\mathrm{K}}(\mathrm{1}-h)({\mathit{\delta }}_{\text{xyl}}+\mathrm{1})+h_{\text{atm}}({\mathit{\delta }}_{\text{atm}}+\mathrm{1})\right]\\ & -\mathrm{1}.\end{array}\end{array}$$

By (i) noting ${\mathit{\alpha }}_{\text{eq}}={\mathit{\epsilon }}_{\text{eq}}+\mathrm{1}$ and ${\mathit{\alpha }}_{\mathrm{K}}={\mathit{\epsilon }}_{\mathrm{K}}+\mathrm{1}$, where ε_eq and ε_K (–) are the equilibrium and kinetic fractionations, respectively, and (ii) dropping terms with products ε_eq⋅ε_K, we obtain the following expression of the difference in isotopic composition between leaf and source waters at ISS:

$$\begin{array}{}\text{(27')}& {\mathit{\delta }}_{\text{L\_ISS}}-{\mathit{\delta }}_{\text{xyl}}={\mathit{\epsilon }}_{\text{eq}}+{\mathit{\epsilon }}_{\mathrm{K}}+h_{\text{atm}}({\mathit{\delta }}_{\text{atm}}-{\mathit{\delta }}_{\text{xyl}}-{\mathit{\epsilon }}_{\mathrm{K}}.)\end{array}$$

We note that Eq. (27') is the inversion of the Craig and Gordon equation at ISS, i.e. when δ_T=δ_xyl. Finally, δ_E is computed with the NSS Craig and Gordon equation, i.e. Eq. (18'). Equation (26) states that, at a permanent state for transpiration, the degree of attainment of ISS conditions in the leaf is a function of time, leaf internal dynamics (τ_L), and (isotopic) aerodynamic boundary conditions. The formula of Dongmann et al. (1974) requires two additional parameters as compared to the more “straightforward” application of the Craig and Gordon equation, namely leaf transpiration (^jF_T) and volume (V_L), both labour-intensive to obtain and associated with high uncertainties.

Both case scenarios (i) and (ii) make the assumption that leaf water is a well-mixed reservoir, in other words that only convective transport of the water isotopologues occurs, leading to δ_L=δ_Lts, where δ_Lts is the isotopic composition of water at the leaf transpiration sites. However, a number of studies reported strong isotopic variations within the leaf water pool (i.e. among different compartments such as leaf veins, cell walls, and symplastic water; see e.g. Yakir et al., 1989; Wang et al., 1998, 1994; Bariac et al., 1994), which can be related to hydraulic separation of water pools and diffusive transport from the transpiration sites towards the petiole of the leaf. Another explanation may be found in the heterogeneity in opening of the leaf stomata (Farquhar et al., 2007). More specifically, δ_Lts should be significantly higher than the bulk leaf water isotopic composition value, δ_L, which leads to an underestimation of δ_E by the direct application of the Craig and Gordon equation. Walker et al. (1989), Walker and Brunel (1990), and Flanagan et al. (1991) considered in a first approach two distinct water pools in the leaf, one in isotopic equilibrium with water vapour in the stomatal cavity (of isotopic composition δ_{Lts_ISS}) and one isotopically undistinguishable from xylem water (of isotopic composition δ_xyl) in respective proportions p and (1−p). In these three studies, an analogous expression to Eq. (27') is used where p is accounted for:

$$\begin{array}{}\text{(27'')}& {\mathit{\delta }}_{\text{Lts\_ISS}}-{\mathit{\delta }}_{\text{xyl}}={\displaystyle \frac{{\mathit{\delta }}_{{\mathrm{L}}_{\text{ISS}}}-{\mathit{\delta }}_{\text{xyl}}}{p}}={\mathit{\epsilon }}_{\text{eq}}+{\mathit{\epsilon }}_{\mathrm{K}}+h_{\text{atm}}({\mathit{\delta }}_{\text{atm}}-{\mathit{\delta }}_{\text{xyl}}-{\mathit{\epsilon }}_{\mathrm{K}}).\end{array}$$

They suggested that there was a midday maximum for T density rate from the corresponding minimum value for p. Cernusak et al. (2002) and Farquhar and Cernusak (2005) proposed a similar equation to that of Dongmann et al. (1974) for the evaporative isotopic enrichment in leaves in NSS conditions but without considering the leaf water volume per unit area constant in time. Equation (25) becomes in their case

$$\begin{array}{}\text{(25')}& \mathrm{d}\left(V_{\mathrm{L}}{\mathit{\delta }}_{\mathrm{L}}\right)={}^{\mathrm{j}}{F}_{\mathrm{T}}\left({\mathit{\delta }}_{\text{xyl}}\right(t)-{\mathit{\delta }}_{\mathrm{T}}(t\left)\right)\mathrm{d}t.\end{array}$$

By replacing δ_xyl and δ_E in the right hand-term of Eq. (25') with the ISS and NSS Craig and Gordon equation forms, respectively, the authors give an expression relating the rate of change of δ_L with the difference between δ_{Lts_ISS}and δ_Lts:

$$\begin{array}{}\text{(28)}& {\displaystyle \frac{\mathrm{d}\left(V_{\mathrm{L}}{\mathit{\delta }}_{\mathrm{L}}\right)}{\mathrm{d}t}}={\displaystyle \frac{{}^{\mathrm{j}}{\mathit{\chi }}_{\text{int}}}{{}^{\mathrm{j}}r\cdot {\mathit{\alpha }}_{\mathrm{K}}{\mathit{\alpha }}_{\text{eq}}}}({\mathit{\delta }}_{\text{Lts\_ISS}}-{\mathit{\delta }}_{\text{Lts}}),\end{array}$$

where ^jχ_int and ^jr are the water vapour mixing ratio in the intercellular space and, as in Sect. 3.2, the resistance to vapour flow of the ${}^{\mathrm{1}}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{16}}\mathrm{O}$ isotopologue in air, respectively. It is therefore possible, by fitting the time course of the bulk leaf water isotopic composition δ_L to deduce δ_Lts, on the basis of which δ_E is finally determined using Eq. (18') (Yepez et al., 2005). α_eq is, as in Sect. 3.2, calculated following the closed-form equations of e.g. Horita and Wesolowski (1994) (Eq. 17). As for α_K, its expression is adapted to include the series of flow resistances of water vapour isotopologues inside the stomatal cavity or through the stomatal opening (ⁱr_sto and ^jr_sto, T L⁻¹) and in the leaf boundary layer (ⁱr_bdl and ^jr_bdl, T L⁻¹) (Jarvis, 1976; Stewart, 1988). Farquhar et al. (1989) (and see also Cernusak et al., 2005; Farquhar et al., 2007) considered that molecular diffusion drives the transport of the different water vapour isotopologues in the first case and that turbulence prevails in the second, leading to n exponent values of 1 and $\mathrm{1}/\mathrm{2}$, respectively (Dongmann et al., 1974; Eq. 19'). In this framework, α_K is decomposed as

$$\begin{array}{}\text{(29)}& \begin{array}{rl}{\mathit{\alpha }}_{\mathrm{K}}={\displaystyle \frac{{}^{\mathrm{i}}r}{{}^{\mathrm{j}}r}}& ={\displaystyle \frac{{}^{\mathrm{i}}{r}_{\text{sto}}+{}^{\mathrm{i}}{r}_{\text{bdl}}}{{}^{\mathrm{j}}{r}_{\text{sto}}+{}^{\mathrm{j}}{r}_{\text{bdl}}}}\\ & ={\displaystyle \frac{{\left(\frac{{}^{\mathrm{j}}D}{{}^{\mathrm{i}}D}\right)}^{\mathrm{1}}\cdot {}^{\mathrm{j}}{r}_{\text{sto}}+{\left(\frac{{}^{\mathrm{j}}D}{{}^{\mathrm{i}}D}\right)}^{\mathrm{1}/\mathrm{2}}\cdot {}^{\mathrm{j}}{r}_{\text{bdl}}}{{}^{\mathrm{j}}{r}_{\text{sto}}+{}^{\mathrm{j}}{r}_{\text{bdl}}}}\end{array}\end{array}$$

Cuntz et al. (2007) proposed a general iterative solution of the Dongmann et al. (1974) formulation revisited by Cernusak et al. (2002) (Eq. 28) under various scenarios depending on considerations regarding leaf water reservoir isotopic homogeneity (δ_L=δ_Lts or δ_L≠δ_Lts) and volume ($\mathrm{d}{V}_{\mathrm{L}}/\mathrm{d}t=\mathrm{0}$ or $\mathrm{d}{V}_{\mathrm{L}}/\mathrm{d}t\ne \mathrm{0}$). Dubbert et al. (2013) applied their solution in the case of an isotopically well-mixed leaf water pool transpiring at constant volume and expressed the incremental change in δ_L from time step t to t+dt as

$$\begin{array}{}\text{(30)}& \begin{array}{rl}{\mathit{\delta }}_{\mathrm{L}}(t+\mathrm{d}t)=& {\mathit{\delta }}_{\text{L\_ISS}}+\left({\mathit{\delta }}_{\mathrm{L}}\right(t)-{\mathit{\delta }}_{\text{L\_ISS}})\\ & \cdot \exp \left(-{\displaystyle \frac{g_{\mathrm{s}}{}^{\mathrm{j}}{\mathit{\chi }}_{\text{int}}}{{\mathit{\alpha }}_{\mathrm{K}}{\mathit{\alpha }}_{\text{eq}}{V}_{\mathrm{L}}}}\mathrm{d}t\right),\end{array}\end{array}$$

where g_s (L T⁻¹) is the total stomatal conductance.

3.3.2 Progress and challenges

The isotopic composition of T may be derived under NSS conditions from plant chamber measurements following Eq. (7) (Sect. 3.1), either at the leaf level or at the branch level. While most studies developed and operated custom-made chambers, only a few (e.g. Wang et al., 2010) used commercially available leaf chambers (e.g. LI-6400, LI-COR, Nebraska, USA). Chamber measurements have several disadvantages as discussed in Sect. 3.1 but are essential for monitoring δ_E directly without relying on additional modelling steps using either δ_xyl or δ_L. The important two features of the chamber-based method are that it does not require the assumption of ISS and that it allows for repeated (i.e. non-destructive) measurements on the same leaf or ensemble of leaves during the course of the day. On the other hand, the determination of δ_xyl or δ_L, which is largely based on destructive sampling and water recovery with e.g. cryogenic vacuum extraction, is also associated with uncertainty (e.g. Orlowski et al., 2016a, b; Millar et al., 2018). The choice of an appropriate method for sampling xylem water is also crucial for a correct determination of δ_E. For example, herbaceous, grass, or crop species do not have suberized stems; thus destructive sampling would have to rely on leaf water sampling or sampling the plant culm belowground, which is highly destructive and may not be possible on plots of common size. Recently, Chen et al. (2020) documented during a series of laboratory controlled experiments that the apparent offset measured between the hydrogen isotopic composition in sap xylem and source water of different mangrove plant species was the result of artefacts during the vacuum extraction process rather than due to isotopic fractionation during water uptake. This could be a reason for hydrogen isotopic offsets reported elsewhere in the literature (e.g. Ellsworth and Williams, 2007; Barbeta et al., 2019). If applicable to other species, the results of Chen et al. (2020) would suggest caution in determining $T/ET$ values based on the determination of δ²H_T directly from δ²H_xyl at ISS (Eq. 24b) or considering NSS conditions using e.g. the Craig and Gordon (1965) equation (Eq. 18').

A novel type of non-destructive method, first published by Volkmann et al. (2016) and recently by Marshall et al. (2020), could enable monitoring δ_E of trees at an equivalent temporal resolution and even greater temporal coverage than with leaf- or plant-scale chamber systems. In the former study, several 10 mm outer diameter gas probes designed after Volkmann and Weiler (2014) (see Sect. 3.2) were inserted into pre-drilled holes in the trunk sapwood of two individuals of Acer campestre L. The probes were positioned at breast height in various azimuths. By assuming isotopic equilibrium between the water vapour sampled by the probe and flushed to the laser spectrometer and the xylem (liquid) water, the authors computed δ_xyl values from the temperature-dependent relationships given by e.g. Horita and Wesolowski (1994). For comparison, tree sapwood was destructively sampled, and its water isotopic composition was measured with IRMS after cryogenic vacuum extraction. A good agreement was found between online measurements and offline analysis of xylem water hydrogen isotopic composition. The inter-method bias regarding the determination of xylem water δ¹⁸O was thought to be due to spectral interferences during online analysis with the laser spectrometer. The experimental natural conditions did not allow the authors to conclude if differences in δ_xyl among the different gas probes reflected actual diurnal variations in root water uptake or preferential connection between xylem vessels and specific parts of the root system that were not affected by the labelling pulse. The authors underline the difficulty with their experimental design to precisely measure the temperature of equilibration in the gas probe (needed for converting sample water vapour to xylem water isotopic composition), due to the high lateral temperature gradient and its daily course. Marshall et al. (2020) tested a cruder way (which they entitled the Borehole Equilibration) to sample water vapour originating from xylem water of two pine tree species (Pinus sylvestris L. and Pinus pinea L.) under semi-controlled conditions. Contrary to Volkmann et al. (2016), (i) the authors did not use a gas probe but simply connected a hole drilled horizontally through the trunk to a laser spectrometer with gas sampling lines. Furthermore, (ii) the experiments were performed in hydroponic water solutions to enable a quasi-instantaneous change of the isotopic composition of the water source, thereby setting defined lower isotopic boundary conditions for further modelling efforts. To test the practicability of the method, the experimental results were confronted with a “Dongmann-like” NSS formulation of the isotopic composition of the water vapour stream, in which the geometry and its consequence on the diffusion from the borehole surface and on the establishment of laminar-flow transport were explicitly accounted for. With their model, the authors tested whether the sampled water vapour was in isotopic equilibrium with xylem water or was the product of evaporation from it. It was shown that the prevalence of a full isotopic equilibrium was a reasonable assumption and that the flow-through time (i.e. borehole volume divided by the flow rate) was 20 times greater than the time needed for diffusion of water vapour originating from the xylem vessels into the laminar-flow region in the middle of the borehole section. Both methods present a drastic advancement in isotopic analysis of xylem water and have great potential in the context of ET partitioning of forest ecosystems, on the pivotal condition that the steady-state assumption (δ_xyl=δ_T) applies during periods of measurements. The long-term applicability of the method, i.e. the ability of the investigated tree species to withstand the invasive and destructive installation of the probe, still needs to be proven at this point.

While the coupling between gas exchange chambers and laser spectrometers has the advantage of directly measuring δ_E, the aforementioned destructive sampling method and in situ monitoring technique quantify δ_L or δ_xyl, therefore potentially requiring a modelling step to obtain δ_E. A number of studies (e.g. Zhou et al., 2018; Wei et al., 2015; Aouade et al., 2016; Volkmann et al., 2016) assume ISS and hence treat δ_xyl as equal to δ_E. Although this assumption is probably justified for a daily integration, there is growing evidence that plants reach ISS only briefly in the course of a day, especially when environmental conditions change rapidly (Simonin et al., 2013; Dubbert et al., 2014b, 2017). Thus the analysis is greatly simplified by daily integration, if that is sufficient for the study objectives. Moreover, the leaf water turnover time, which can effectively be described by stomatal conductance (g_s), vapour pressure deficit, and leaf water volume, is species-specific and ranges from several minutes to several hours (Song et al., 2015). As the leaf water turnover time describes the necessary time for leaf water to reach ISS (see exponent terms in Eqs. 26 and 30), ISS can either be observed for large parts of the day (e.g. in many herbaceous species) or not at all (e.g. in plant species strongly controlling their g_s; see Dubbert et al., 2017, and Dubbert and Werner, 2019, for an overview). Therefore, the validity of assuming ISS for the purpose of ET partitioning will largely depend on the desired temporal scale: considering NSS is definitely necessary at subdiurnal to diurnal scales but unimportant at larger timescales. In case NSS is likely to occur, δ_E can be modelled using a “Dongmann version” of the Craig and Gordon equation, as shown in Sect. 3.3.1 (Dongmann et al., 1974). However, this complicates the partitioning approach considerably in comparison to direct chamber measurements of δ_E, as a large number of additional observations are necessary. In particular, g_s and the canopy temperature are important input parameters. Therefore, the use of chamber measurements is highly recommended in any case.

4 Summary and possible ways forward

The isotopic methodology for partitioning ET relies on a number of possible combinations of different techniques, which differ in numerous aspects. While some of them are based on destructive sampling and water recovery using one of the aforementioned methods (e.g. cryogenic vacuum extraction or direct liquid–vapour equilibration; see Sect. 3.2) and a posteriori analysis in the laboratory (e.g. for determination of δ_E using the Craig and Gordon equation), other methods are non-destructive, provide online measurements, and do not include a strong modelling component (e.g. determination of δ_E with plant chambers with one of two mass balance techniques). Destructive approaches require neither the handling of soil, plant, or soil and plant chambers nor the deployment of a laser spectrometer along with its conditioning system in the field. They should also allow for capturing the inherent spatial variability with repeated sampling (however, at the cost of long hours spent in sample preparation and water extraction). Non-destructive methods, such as chambers, may on the other hand provide environmental conditions for the enclosed plant that are not representative of ambient conditions.

Up to now only indirect methods, e.g. based on Scanlon and Kustas (2010), might have been able to provide continuous and subdaily estimates of $T/ET$. Some methods, such as the Keeling plot technique, can provide long-term continuous estimates of δ_ET once a meteorological mast is installed in the field. It is, on the other hand, not advisable to enclose a plant in a chamber over longer time periods. Within the realm of destructive techniques, the user may assume ISS or test its existence when determining the isotopic composition of T. The techniques with which δ_ET is estimated generally differ in terms of spatial significance as compared to those for determining δ_E and δ_E. Estimates of δ_ET obtained either with the EC, Keeling plot, or flux gradient technique are thought to be representative at the field scale (e.g. as represented by the EC footprint). Note that this is also a problem encountered in (non-isotopic) instrumental approaches for partitioning ET, including EC, micro-lysimeters, and soil chambers (Kool et al., 2014). To account for these discrepancies in spatial representativeness, several micro-lysimeters and (if possible automated) chambers are deployed on site, e.g. within the framework of global networks (e.g. FLUXNET; Law et al., 2002). On the contrary, there has been no consensus to date on a common methodological ground for partitioning ET in the field on the basis of water stable isotopic measurements, depending on the type of land cover and use (agricultural, grassland, or forest ecosystems).

Figure 6 One example of application of the Keeling (1958) plot technique to synthetic data that would be collected with a field-deployable lift at high vertical resolution (0.1 m) (for implications on measurement frequency, which also needs to be high (≥ 5 Hz), see Sects. 3.1 and 4 of Ney and Graf, 2018). The oxygen isotopic compositions of evapotranspiration, transpiration, and evaporation are estimated by the values of the y intercepts of the linear regressions between the isotopic composition of the atmospheric water vapour (δ_atm) and the inverse of the water vapour mixing ratio (χ_atm) in three non-overlapping regions, i.e. (i) the “free atmosphere” (indicated by the blue symbols), (ii) the region spreading from the canopy height to the height of the local maximum in δ_atm (green symbols), and (iii) the region delimited by the δ_atm local maximum height and the ground level (brown symbols). Also shown are the 95 % confidence interval envelopes of the linear regressions (dashed lines) as well as error bars (1 standard error) of the y intercepts.

It is the authors' belief that non-destructive and online methods integrated into automated sampling platforms and part of long-term (e.g. multi-year) water flux observatories should be preferred over destructive and discontinuous assessments of $T/ET$ values. In this (ideal) framework, we propose the following:

• i.

  The seminal effort in applying the EC technique by Griffis et al. (2010) should be continued to provide half-hourly and continuous ecosystem-scale δ_ET estimates. The δ_ET estimates obtained with the EC technique should be corroborated or confronted with the Keeling plot and the flux gradient approaches to identify possible scale-dependent disparities in surface isotopic signals as in Good et al. (2012).

• ii.

  δ_E should be monitored by installing gas-permeable membranes or tubing (see Sect. 3.2) in the upper layers of the soil, depending on site-specific knowledge regarding the receding of the EF. While the gas probes of Volkmann and Weiler (2014) and Gaj et al. (2016) are better-suited for insertion at different locations in a soil profile, the membrane tubing used by Rothfuss et al. (2013), Oerter and Bowen (2019), and Kübert et al. (2020) allows for covering more ground surface by using a customized length of tubing. This should help to increase the representativeness of the δ_E value estimated from the soil water vapour isotopic composition and the use of the Craig and Gordon equation. When using the model of Craig and Gordon (1965), authors should systematically perform sensitivity analyses of

  • a.

    the depth of the EF and its water isotopic composition and

  • b.

    the value of the kinetic fractionation factor, α_K.

  These analyses will provide insights into the uncertainty of $T/ET$, in addition to the uncertainty originating from the solution of the two end-member equation (Eq. 1) (Rothfuss et al., 2010). This is, however, under-investigated according to our literature review. The α_K value may be derived in a dual-isotope space using the formulation of Gat (2000) rather than based on unclear assumptions regarding the type of transport (molecular diffusion, laminar, or turbulent transport) controlling the flow of water stable isotopologues (see Sect. 3.2). As a side note (and without a proof of concept for this), the δ_E value may be directly determined in the case of a well-developed dry surface on the basis of non-destructive measurements of the soil water vapour isotopic composition (${\mathit{\delta }}_{\text{soil}}^{\text{vap}}$) at two depths (z₁ and z₂) located between the EF and the soil surface. For this, the Craig and Gordon equation may be used without the need to locate the soil EF nor to assume liquid–vapour equilibrium:

  $$\begin{array}{}\text{(18'')}& {\mathit{\delta }}_E={\displaystyle \frac{{\mathit{\delta }}_{\text{soil}}^{\text{vap}}\left(z_{\mathrm{1}}\right)\cdot h\left(z_{\mathrm{1}}\right)-{\mathit{\delta }}_{\text{soil}}^{\text{vap}}\left(z_{\mathrm{2}}\right)\cdot h\left(z_{\mathrm{2}}\right)}{\frac{{}^{\mathrm{j}}D}{{}^{\mathrm{i}}D}\left(h\right(z_{\mathrm{1}})-h(z_{\mathrm{2}}\left)\right)}}-\mathrm{1}.\end{array}$$
• iii.

  Several transparent flushed plant- or leaf-size chambers should be operated at the study site to characterize the in situ natural lateral heterogeneity of δ_E, due to differences in root water uptake, plant physiological state, and lateral heterogeneity in soil water isotopic composition profiles. This would be a prerequisite for any upscaling attempt of δ_E values. Developments should be made towards designing chambers able to mimic the dynamics of ambient air temperature, relative humidity, and wind turbulence to avoid biases in δ_E estimation. This could be done by cooling of the inlet air to avoid overheating of the air inside the chamber and with an adaptive active ventilation system. In situations where parts of the field are bare, e.g. between crop rows, soil chambers should be installed as well to evaluate differences in δ_E between areas covered or not covered with vegetation.

• iv.

  The methods for monitoring of δ_xyl and its potential use in determining δ_E (that is by assuming ISS conditions) have been tested and validated with tree species exclusively. The same principle is yet to be minimized and applied to crops able to survive the installation and carry the instrumentation, such as a well-developed maize plant.

Lastly, the lift system principle, as operated by Noone et al. (2013), Mayer et al. (2009), and recently for agricultural crops Ney and Graf (2018) has the potential to provide half-hourly concomitant values of δ_ET, δ_E, and δ_E in the field. The principle is illustrated in Fig. 6, further developing that of Yepez et al. (2003). The Keeling plot technique is applied to data collected at high vertical resolution (ultimately implying high-frequency data acquisition of the analyser, typically equal to or greater than 5 Hz; see Ney and Graf, 2018) in three distinct atmospheric regions, i.e. (i) the region spreading from the fully turbulent atmosphere to the canopy height, (ii) the region comprised between canopy height (here fixed at 1.25 m) and the local maximum in δ_atm, and (iii) the region delimited by the δ_atm local maximum and the ground level (Fig. 6a). The y intercepts of the three Keeling plots give the concomitant values of the isotopic compositions of ET (Fig. 6b), T (Fig. 5c), and E (Fig. 6d). In this synthetic experiment, which cannot be construed as a proof of concept, δ¹⁸O_ET, δ¹⁸O_T, and δ¹⁸O_E are equal to −4.7 (± 1.5), −0.7 (± 1.4), and −18.5 (± 0.4) ‰, respectively, corresponding to a $T/ET$ value of 77 (± 10) %.

Importantly, and to conclude this summary, the general isotopic partitioning approach (i.e. Eq. 1) as well as the ensemble of methods and their possible improvements will not be applied broadly until they are able to deal with canopy interception. Further research is therefore needed to (i) determine the water volumes collected by the vegetation following rain events, fog deposition, or dew condensation and (ii) investigate the isotopic effects during re-evaporation of the intercepted water. This should be useful for constraining a generalized partitioning equation including a third member, namely the stable isotopic composition of interception.

5 Conclusions

Water stable isotopes are often described in the present literature compilation as “powerful” (or “insightful”) tools for separating evaporation and transpiration fluxes. However, the number of ET partitioning studies, which the authors listed here, remains low when compared to the number of publications utilizing water stable isotopes for e.g. determining plant water use strategies (30 vs. 158 over the period 1990–2016; see Rothfuss and Javaux, 2017). The apparent contrast between the announced potential and the number of study cases is explained partly by both the complexity and multifaceted character of the isotopic methodology. Unfortunately, despite great efforts of the researchers, the spatial representativeness and temporal extent of the obtained $T/ET$ data series are usually not well comparable with those of other non-isotopic methods (see Fig. 2g).

The authors believe that while ultimately increasing the complexity in terms of modus operandi, novel non-destructive monitoring methods are key to providing long-term $T/ET$ data at the plot to the field scale and to upscaling local process understanding to address large-scale ecohydrological issues in a changing climate.

Appendix A

Table A1 Overview of the partitioning studies found with the ISI Web of Science search engine (http://www.webofknowledge.com/, last access: 15 February 2021) on basis of search terms ((“evapotranspiration” or “transpiration” or “evaporation”) and partition* and isotop*). “CG65”, “Kp58”, and “f-g” refer, respectively, to the Craig and Gordon (1965) equation for determination of δ_E and δ_E, the Keeling (1958) plot, and flux gradient techniques for determination of δ_ET and δ_E. “Chamber(InOut)” and “Chamber(Kp58)” refer to gas-exchange-chamber-based measurements for determination of δ_ET, δ_E, and δ_E either by comparing the chamber inlet and outlet gas properties or by applying the Keeling plot technique, respectively. The reader is referred to the list of symbols and abbreviations in Appendix B for the definitions of the other symbols and abbreviations.

Appendix B: List of symbols and abbreviations used in the main document and Table A1

Symbol or abbreviation	Description	Dimension or unit
Catm, Cbg, CET	Atmospheric and background water vapour concentration, rise in atmospheric water vapour concentration due to evapotranspiration flux	M L−3
E,ELys,Epot	Soil evaporation rate, soil evaporation (micro-lysimeter), potential evaporation	L3 T−1
EC	Eddy covariance
ET	Evapotranspiration rate	L3 T−1
f	Measurement frequency	T−1
FET, jFET, iFE, FE, jFE, jFT, iFT, jFxyl, iFxyl	Evapotranspiration water vapour flux density rate, evapotranspiration, evaporation, transpiration, and xylem water flux density rates of the rare (i) and abundant (j) isotopologue	${\mathrm{L}}^{\mathrm{3}}{\mathrm{L}}^{-\mathrm{2}}{\mathrm{T}}^{-\mathrm{1}}$
gs	Leaf stomatal conductance	$\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$
GPP	Gross primary production	$\mathrm{M}{\mathrm{L}}^{-\mathrm{2}}{\mathrm{T}}^{-\mathrm{1}}$
hatm,hEF	Atmospheric relative humidity and soil pore space relative humidity at the evaporating front	–
K	Eddy diffusivity of water vapour	L2 T−1
LAI	Leaf area index	L2 L−2
LET	Latent heat flux of evapotranspiration	M T−3
Matm, Mw	Molecular weight of dry air and water	M L−3
n	Adimensional factor accounting for flow conditions above the liquid water–water vapour equilibrium layer	–
NEE	Net ecosystem exchange	$\mathrm{M}{\mathrm{L}}^{-\mathrm{2}}{\mathrm{T}}^{-\mathrm{1}}$
P	Precipitation amount	L3 L−2
p	Proportion of leaf water in isotopic equilibrium with water vapour in the stomatal cavity	–
PPFD	Photosynthetic photon flux density	$\mathrm{\mu }\mathrm{mol}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{2}}$
Qs	Sensible heat flux	M T−3
R	Universal gas constant	$\mathrm{M}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{T}}^{-\mathrm{3}}$
Rstd, REF, Rsat	Isotopic ratio of the Vienna Standard Mean Ocean Water (V-SMOW), soil water at the evaporating front, and saturated water vapour	–
Rn, Rs, Rd	Net and solar radiation and radiation flux density	M T−3
Rgas	Universal gas constant	$\mathrm{M}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{T}}^{-\mathrm{3}}$
jr, ir jD, iD	Bulk resistances to vapour transport of the rare (i) and abundant (j) isotopologuesMolecular diffusivities of the rare (i) and abundant (j) water vapour isotopologues
S	Sap-flux density	$\mathrm{M}{\mathrm{L}}^{-\mathrm{2}}{\mathrm{S}}^{-\mathrm{1}}$
T	Transpiration rate	L3 T−1
Tatm, Tsoil, TEF, TL, Tcan	Temperature of the atmosphere, soil, soil at the evaporating front, leaf surface, and canopy atmosphere	∘C
$T/ET$	Transpiration fraction	–
uin, uout	Flow rate measured at the inlet and outlet of a gas exchange chamber	L3 T−1
v (vd)	Wind speed (wind direction)	L T−1
VFC	Vegetation fractional coverage	L2 L−2
VPD	Vapour pressure deficit	P
z, zatm, zEF	Vertical coordinate, atmospheric height, and soil evaporating front depth	M
αeq αK	Equilibrium isotopic fractionation factorKinetic-isotopic-fractionation factor	–
δatm, δsoil, ${\mathit{\delta }}_{\text{soil}}^{\text{vap}}$, δL, δxyl, δprec, δroot, δirr, δpond, δin, δout	Isotopic composition of the atmospheric water vapour, soil water, soil water vapour, leaf water, xylem water, precipitation, root water, irrigation water, pond water, and water vapour measured at the inlet and outlet of a gas exchange chamber	–
εeq εK	Equilibrium isotopic fractionationKinetic isotopic fractionation	–
φ	Isotope analyser inlet flow rate	L3 T−1
ρatm, ρw	Volumetric mass of dry air and water	M L−3
${\mathit{\theta }}_{\text{soil}},{\mathit{\theta }}_{\text{surf}},{\mathit{\theta }}_{\text{res}},{\mathit{\theta }}_{\text{sat}},{\mathit{\theta }}_{\mathrm{L}}$	Soil, soil surface, soil residual, and soil saturated water content, leaf water content	L3 L−3
τL	Leaf water turnover time	T
χatm, χbg, jχatm, iχatm, ${}^{\mathrm{j}}{{\mathit{\chi }}_{\text{atm}}}^{\text{sat}}$, ${}^{\mathrm{i}}{{\mathit{\chi }}_{\text{atm}}}^{\text{sat}}$, χET, χin, χout	Atmospheric and background water vapour mixing ratio, water vapour mixing ratio of the rare (i) and abundant (j) isotope, saturated water vapour mixing ratio of the rare (i) and abundant (j) isotope, rise in atmospheric water vapour mixing ratio due to evapotranspiration flux, water vapour mixing ratio measured at the inlet and outlet of a gas exchange chamber	L3 L−3
ψEF	Soil water matric potential at the evaporating front	$\mathrm{M}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{T}}^{-\mathrm{2}}$