Reviews and syntheses: Turning the challenges of partitioning ecosystem evaporation and transpiration into opportunities

Evaporation (E) and transpiration (T) respond differently to ongoing changes in climate, atmospheric composition, and land use. It is difficult to partition ecosystem-scale evapotranspiration (ET) measurements into E and T, which makes it difficult to validate satellite data and land surface models. Here, we review current progress in partitioning E and T and provide a prospectus for how to improve theory and observations going forward. Recent advancements in analytical techniques create new opportunities for partitioning E and T at the ecosystem scale, but their assumptions have yet to be fully tested. For example, many approaches to partition E and T rely on the notion that plant canopy conductance and ecosystem water use efficiency exhibit optimal responses to atmospheric vapor pressure deficit (D). We use observations from 240 eddy covariance flux towers to demonstrate that optimal ecosystem response to D is a reasonable assumption, in agreement with recent studies, but more analysis is necessary to determine the conditions for which this assumption holds. Another critical assumption for many partitioning approaches is that ET can be approximated as T during ideal transpiring conditions, which has been challenged by observational studies. We demonstrate that T can exceed 95 % of ET from certain ecosystems, but other ecosystems do not appear to reach this value, which suggests that this assumption is ecosystem-dependent with implications for partitioning. It is important to further improve approaches for partitioning E and T, yet few multi-method comparisons have been undertaken to date. Advances in our understanding of carbon–water coupling at the stomatal, leaf, and canopy level open new perspectives on how to quantify T via its strong coupling with photosynthesis. Photosynthesis can be constrained at the ecosystem and global scales with emerging data sources including solar-induced fluorescence, carbonyl sulfide flux measurements, thermography, and more. Such comparisons would improve our mechanistic understanding of ecosystem water fluxes and provide the observations necessary to validate remote sensing algorithms and land surface models to understand the changing global water cycle.

1 Introduction

Some 70 000 km³ of water leaves terrestrial ecosystems and enters the atmosphere through evapotranspiration (ET) every year (Jung et al., 2019; Oki and Kanae, 2006). Despite its importance, we are unsure whether global ET has been increasing over time (Brutsaert, 2013, 2017; Brutsaert and Parlange, 1998; Zeng et al., 2018; Zhang et al., 2016) such that the water cycle is accelerating (Ohmura and Wild, 2002) or decreasing and causing more river discharge (Gedney et al., 2006; Labat et al., 2004; Probst and Tardy, 1987). Global ET volumes from reanalyses, upscaled estimates, and land surface model (LSM) outputs disagree (Mueller et al., 2013) by up to 50 % (Mao et al., 2015; Vinukollu et al., 2011). LSMs also struggle to simulate the magnitude and/or seasonality of ET at the ecosystem scale (Fig. 1), suggesting fundamental gaps in our understanding of the terrestrial water cycle. These issues need to be resolved to effectively manage water resources as climate continues to change (Dolman et al., 2014; Fisher et al., 2017).

Figure 1 The mean monthly latent heat flux (λE) – the energy used for evapotranspiration – from eddy covariance measurements from four research sites (“MEASURED”) and 13 ecosystem models from the North American Carbon Program Site-Level Interim Synthesis (Schwalm et al., 2010). Sites: CA-Ca1 (Schwalm et al., 2007), CA-Obs (Griffis et al., 2003; Jarvis et al., 1997), US-Ho1 (Hollinger et al., 1999), US-Me2 (Thomas et al., 2009). Models: BEPS (Liu et al., 1999), CAN-IBIS (Williamson et al., 2008), CNCLASS (Arain et al., 2006), ECOSYS (Grant et al., 2005), ED2 (Medvigy et al., 2009), ISAM (Jain and Yang, 2005), ISOLSM (Riley et al., 2002), LOTEC (Hanson et al., 2004), ORCHIDEE (Krinner et al., 2005), SIB (Baker et al., 2008), SIBCASA (Schaefer et al., 2009), SSIB2 (Zhan et al., 2003), TECO (Weng and Luo, 2008). Data are available from Ricciuto et al. (2013).

Along with technological and data limitations, we argue that a fundamental challenge in modeling ET at the global scale is difficulty measuring transpiration (T) through plant stomata and evaporation (E) from non-stomatal surfaces at the ecosystem scale (Fisher et al., 2017; McCabe et al., 2017). LSMs and remote sensing algorithms (see Appendix A) rely on a process-based understanding of E and T to estimate ET, but it is not clear how to guide their improvement without accurate ground-based E and T observations at spatial scales on the order of a few kilometers or less (Talsma et al., 2018) and temporal scales that capture diurnal, seasonal, and interannual variability in water fluxes. Recent statistical ET partitioning approaches (Rigden et al., 2018) are similarly limited by the lack of direct E and T observations for evaluation. Interest in partitioning E and T from ecosystem ET measurements has grown in recent years (Anderson et al., 2017b), and many new measurements and modeling approaches seek to do so but often rely on assumptions that need further testing. We begin with a brief research review that notes recent updates to our theoretical understanding of ET and outlines the challenges in measuring E and T at the ecosystem scale. We then describe current and emerging innovations in partitioning E and T (Table 1) and use observations to challenge some of the assumptions upon which these approaches rely. We finish with an outlook of how carefully designed ecosystem-scale experiments can constrain models of E and T to improve our understanding going forward.

Table 1 A summary of recent approaches for estimating transpiration (T) and/or for partitioning evapotranspiration (ET) into evaporation (E) and T at the ecosystem scale. The reader is referred to Kool et al. (2014) for a comprehensive review of E and T measurement methodologies.

2 Background

2.1 Vegetation plays a central role in evaporation and transpiration partitioning

The ratio of transpiration to evapotranspiration (T ∕ ET) at annual timescales is related to aridity (Good et al., 2017) but appears to be relatively insensitive to annual precipitation (P) (Schlesinger and Jasechko, 2014). T ∕ ET is sensitive to ecosystem characteristics, namely the leaf area index (LAI) (Berkelhammer et al., 2016; Fatichi and Pappas, 2017; Wang et al., 2014; Wei et al., 2015), especially on sub-annual timescales (Li et al., 2019; Scott and Biederman, 2017), noting that LAI is related to P at longer timescales. A higher LAI favors T and E from intercepted water (E_i) at the expense of E from soil (E_soil) such that LAI explains some 43 % of the variability of annual T ∕ ET across global ecosystems (Wang et al., 2014). Upscaling this relationship results in a global estimate of terrestrial annual T ∕ ET of 0.57±0.07 (Wei et al., 2017). Other observational studies suggest that annual T ∕ ET averages nearly 2∕3 globally (0.61±0.15, Schlesinger and Jasechko, 2014; 0.64±0.13, Good et al., 2015; and 0.66±0.13 across some FLUXNET sites, Li et al., 2019). Intercomparison studies agree on the large uncertainty surrounding these estimates, with reported global terrestrial annual T ∕ ET ratios ranging from 0.35 to 0.90 (Coenders-Gerrits et al., 2014; Fatichi and Pappas, 2017; Young-Robertson et al., 2018). Approaches that use stable isotopes tend to produce higher annual T ∕ ET values due to assumptions regarding isotopic fractionation (Jasechko et al., 2013; Sutanto et al., 2014). Some LSM estimates of annual T ∕ ET arrive at larger values on the order of 0.70±0.09 (Fatichi and Pappas, 2017; Paschalis et al., 2018), while other LSMs suggest smaller T ∕ ET; for example, T ∕ ET from the IPCC CMIP5 intercomparison ranges from 0.22 to 0.58 (Wei et al., 2017). Constraining these model results with observations results in an estimate similar to observational studies but with reduced uncertainty: 0.62±0.06 (Lian et al., 2018). A number of recent studies suggest that a major cause of the discrepancies between observations and LSM predictions of T ∕ ET is the treatment of lateral flow in models (Chang et al., 2018). Explicitly adding lateral flow and groundwater dynamics is critical for accurate T estimation (Maxwell and Condon, 2016), and realistic lateral flow results in lower E (Chang et al., 2018; Ji et al., 2017). Simulating sub-grid water partitioning is often of particular importance during drought (Ji et al., 2017; Shrestha et al., 2018), as is a realistic representation of plant water stress parameters (Fang et al., 2017). In addition to challenges in simulating T ∕ ET across space, we also need to measure and model T ∕ ET correctly at the ecosystem scale across all timescales over which it varies from minutes or less to multiple years or more. For this, an understanding of ecosystem water transport and biological responses to micrometeorological forcing is necessary (Badgley et al., 2015).

2.2 Turning theory into practice

Measuring and modeling water fluxes from the surface to the atmosphere at the ecosystem scale across multiple scales in time is a nontrivial challenge. The pools in which water is stored in ecosystems span spatial scales from soil pores to forest canopies. Liquid and gaseous water transport occurs through pathways in the soil, xylem, leaves, and plant surfaces that exhibit nonlinear responses to hydroclimatic forcing, which is itself stochastic (Katul et al., 2007, 2012). These complex dynamics of water storage and transport impact the conductance of water between ecosystems and the atmosphere (Mencuccini et al., 2019; Siqueira et al., 2008), and these conductance terms are central to the Penman–Monteith equation, which combines the thermodynamic, aerodynamic, environmental, and biological variables to which ET (m s⁻¹) responds to represent the mass and energy balance of water flux between the land surface and the atmosphere (Monteith, 1965; Penman, 1948):

$$\begin{array}{}\text{(1)}& \mathrm{ET}={\displaystyle \frac{\mathrm{1}}{\mathit{\rho }\mathit{\lambda }}}{\displaystyle \frac{s\left(R_{\mathrm{n}}-G\right)+{\mathit{\rho }}_{\mathrm{a}}{c}_{p}D{g}_{\mathrm{a}}}{s+\mathit{\gamma }\left(\mathrm{1}+\frac{g_{\mathrm{a}}}{g_{\mathrm{surf}}}\right)}}.\end{array}$$

In the Penman–Monteith equation, λ is the latent heat of vaporization (J kg⁻¹), ρ is the density of water (kg m⁻³), s is the slope of the saturation vapor pressure function (Pa K⁻¹), R_n is the surface net radiation (W m⁻²), G is the ground heat flux (W m⁻²), ρ_a is dry air density (kg m⁻³), c_p is the specific heat capacity of air (J kg⁻¹ K⁻¹), D is the vapor pressure deficit (Pa), γ is the psychrometric constant (Pa K⁻¹), g_a is the conductance of the atmosphere, and g_surf is surface conductance to water vapor flux (both m s⁻¹). g_surf is a spatially upscaled effective parameter that includes canopy conductance from stomatal opening (g_c) associated with T, conductance related to soil evaporation (g_soil) associated with E_soil, and conductance related to plant-intercepted evaporation (g_i) associated with E_i. The combination of E_soil and E_i results in ecosystem-scale E. The biological drivers that alter g_c impact T, but physical drivers impact both E and T. In practice, the Penman–Monteith equation is commonly simplified because of the challenge of correctly simulating all relevant conductances (Maes et al., 2019; Priestley and Taylor, 1972).

The micrometeorological drivers of the Penman–Monteith equation vary within and across plant canopies and landscapes (Jarvis and McNaughton, 1986), as do the turbulent structures that transport water into the atmosphere by which ET can be measured using eddy covariance. Because ET is commonly measured above plant canopies with eddy covariance, micrometeorological variables are commonly measured above plant canopies as well. These measurements do not necessarily reflect micrometeorological conditions at evaporating and transpiring surfaces. For example, characteristic profiles of water vapor concentration in the atmosphere measured above the plant canopy are different from D at the canopy, leaf, and soil levels (De Kauwe et al., 2017; Jarvis and McNaughton, 1986; Lin et al., 2018). Furthermore, the fundamental assumption that D reflects the difference between atmospheric water vapor pressure and saturated conditions within the leaf is challenged by studies demonstrating that leaf vapor pressure need not be saturated (Cernusak et al., 2018). Radiation, temperature, and wind speed also vary throughout plant canopies with consequences for modeling T from the canopy and E from the soil and other ecosystem surfaces. The space–time variability of environmental drivers within plant canopies should therefore ideally be measured or simulated to understand how they impact E and T, and ecosystem modelers must decide if this canopy-resolved detail is important to simulate in diverse ecosystems (Boulet et al., 1999; Medvigy et al., 2009; Polhamus et al., 2013).

Modeling ET at the ecosystem scale is challenging enough before noting that ongoing changes to the Earth system impact all of the biotic and abiotic variables that determine it. The decline in incident radiation across some regions of the world due largely to anthropogenic aerosols (“global dimming”) and subsequent increase since about 1990 (“global brightening”) have changed incident radiation and thus R_n at the land surface (Wild et al., 2005). The observed decrease in wind speed (“global stilling”) (McVicar et al., 2012a, b) is partly due to increases in surface roughness owing to increases in LAI (Vautard et al., 2010) and has decreased g_a, which is a function of wind speed (Campbell and Norman, 1998). Atmospheric heating changes the terms in Eq. (1) that involve temperature, namely R_n (via incident longwave radiation), λ, γ, and s, through the Clausius–Clapeyron relation. A warming climate also increases D in the absence of changes in specific humidity, but specific humidity has increased across many global regions (Willett et al., 2008), resulting in complex spatial and temporal changes in D (Ficklin and Novick, 2017). g_c is controlled by soil moisture availability (Porporato et al., 2004), plant hydrodynamics (Bohrer et al., 2005; Matheny et al., 2014), and environmental variables including D that result in stomatal closure (Oren et al., 1999) (Fig. 2), which is critical for models to accurately simulate (Rogers et al., 2017). This dependency on D is predicted to become increasingly important as global temperatures continue to rise (Novick et al., 2016), but D is also highly coupled to soil moisture (Zhou et al., 2019), and both depend on ET itself through soil–vegetation–atmosphere coupling. Increases in atmospheric CO₂ concentration tend to decrease stomatal conductance at the leaf scale (Field et al., 1995) and have been argued to decrease g_c on a global scale (Gedney et al., 2006). However, elevated CO₂ often favors increases in LAI (e.g., Ellsworth et al., 1996), thus leading to an increase in transpiring area that can support greater g_c. Atmospheric pollutants including ozone also impact g_c with important consequences for vegetation function (Hill et al., 1969; Wittig et al., 2007). Water fluxes from the land surface impact atmospheric boundary layer processes including cloud formation, extreme temperatures, and precipitation (Gerken et al., 2018; Lemordant et al., 2016; Lemordant and Gentine, 2018), which feeds back to land surface fluxes in ways that are inherently nonlinear and difficult to simulate (Ruddell et al., 2013). In addition to these highly nonlinear dynamics of the soil–vegetation–atmosphere system, ongoing land use and land cover changes impact vegetation structure and function with important implications for the water cycle. In brief, we need to correctly simulate how E and T respond to a range of biotic and abiotic variability for predictive understanding. To do so, we need to accurately measure E and T in the first place.

Figure 2 The relationship between above-canopy vapor pressure deficit (D) and evapotranspiration (ET in millimeters per half hour, hh) visualized using kernel density estimation (Botev et al., 2010) for more than 1.5 million half-hourly eddy covariance observations with a solar zenith angle less than 60^∘ from 241 eddy covariance research sites in the La Thuile FLUXNET database that included ecosystem type and soil heat flux measurements described in Stoy et al. (2013).

3 Measuring and estimating evaporation and transpiration

There are multiple established methods to measure ecosystem E and T, including leaf gas exchange, plant-level sap flow, lysimeters, soil, leaf, and canopy chambers, photometers, soil heat pulse methods, and stable and radioisotopic techniques. Ongoing efforts to synthesize measurements of ecosystem water cycle components – for example, SAPFLUXNET (Poyatos et al., 2016) – are a promising approach to build an understanding of different terms of the ecosystem water balance across global ecosystems. Multiple reviews and syntheses of E and T measurements have been written (e.g., Abtew and Assefa, 2012; Anderson et al., 2017b; Blyth and Harding, 2011; Kool et al., 2014; Shuttleworth, 2007; Wang and Dickinson, 2012) and have provided the key insights that ecosystem models use to simulate ecosystem–atmosphere water flux (De Kauwe et al., 2013). Rather than reiterate the findings of these studies, we focus on existing and emerging approaches to partition E and T at the ecosystem scale on the order of tens of meters to kilometers at temporal resolutions on the order of minutes to hours, with a particular emphasis on new observational and methodological techniques. We do so to align ecosystem-scale observations of E and T with satellite-based algorithms that can scale E and T from ecosystem to region to globe (Appendix A).

ET is commonly approximated as the residual of the water balance at the watershed scale in hydrologic studies – especially when the change in water storage can be assumed to be negligible – but can now be measured using eddy covariance at the ecosystem scale (Wilson et al., 2001). Other approaches including scintillometry (Cammalleri et al., 2010; Hemakumara et al., 2003), surface renewal (Snyder et al., 1996), and the Bowen ratio energy balance method provide important complements to eddy covariance techniques for measuring ecosystem-scale ET. Such syntheses follow ongoing efforts to compile ET measured by eddy covariance via FLUXNET and cooperating consortia (Chu et al., 2017), which synthesize half-hourly to hourly eddy covariance flux measurements that have been used to partition ET into E and T with mixed success.

3.1 Partitioning ET using half-hourly eddy covariance observations

An early attempt to partition E and T directly from eddy covariance measurements assumed that ET is comprised solely of E in the absence of canopy photosynthesis (gross primary productivity, GPP) due to the coupled flux of carbon and water through plant stomata (Stoy et al., 2006). It was further assumed that E_soil dominated ET during these times and that E_soil could be modeled by simulating solar radiation attenuation through grass, pine forest, and deciduous forest canopies in the Duke Forest, NC, USA. T was subsequently approximated as the difference between measured ET and the model for E_soil during times when photosynthesis was active. Annual T ∕ ET values from this approach varied from 0.35 to 0.66 in the grass ecosystem (US-Dk1) across a 4-year period and between 0.7 and 0.75 in the pine (US-Dk3) and hardwood (US-Dk2) forests, somewhat higher than global syntheses (Schlesinger and Jasechko, 2014), remote sensing estimates from PT-JPL (see Appendix A) for the Duke pine forest (Fig. 3), and sap-flow-based measurements from the deciduous forest (Oishi et al., 2008). These discrepancies arose in part because E_i was considered negligible but can be considerable (see Sect. 3.6). The model for E_soil could also not be directly validated using measurements from the forest floor alone with available observations.

Figure 3 The Priestley–Taylor Jet Propulsion Lab (PT-JPL) estimate of transpiration (T) in energy flux units compared against T estimated using eddy covariance measurements and models of soil evaporation in a loblolly pine forest for 2001–2005 from Stoy et al. (2006). Measurements were taken at 10:30 Eastern Standard Time (UTC−5 h).

An under-explored approach for partitioning E_soil from ecosystem ET uses concurrent above- and below-canopy eddy covariance measurements in forest and savanna ecosystems (Misson et al., 2007). Subcanopy eddy covariance measurements have proven useful for measuring below-canopy ET, often assumed to be comprised largely of E_soil in ecosystems with poor understory cover (Baldocchi et al., 1997; Baldocchi and Ryu, 2011; Moore et al., 1996; Sulman et al., 2016). However, such measurements are not yet widely adopted for ET partitioning studies due to a limited understanding of their performance (Perez-Priego et al., 2017); most work to date has used below-canopy eddy covariance to partition canopy GPP and soil respiration (Misson et al., 2007). Several recent studies demonstrated the additional value of concurrent below-canopy measurements for quantifying the coupling and decoupling of below- and above-canopy airspace to accurately apply the eddy covariance technique in forested ecosystems (Jocher et al., 2017, 2018; Paul-Limoges et al., 2017; Thomas et al., 2013), arguing that below-canopy eddy covariance measurements should be more widely adopted. Other eddy-covariance-based partitioning methods take a different approach and use the relationship between T and GPP to partition ecosystem-scale E and T.

Scott and Biederman (2017) assumed that T is linearly related to GPP at monthly timescales over many years such that

$$\begin{array}{}\text{(2)}& T=m_{\mathrm{WUE}}{r}_{\mathrm{GPP}},\end{array}$$

where m_WUE is the inverse of the marginal water use efficiency (the change, Δ, in ET per change in GPP: ΔET∕ΔGPP), and r is the ratio between the inverse of the transpirational water use efficiency (ΔT∕ΔGPP) and the marginal ecosystem water use efficiency, which is assumed to be unity. It follows that the intercept E^′ of the relationship $\mathrm{ET}=m_{\mathrm{GPP}}+E^{\prime }$ is an estimate of average monthly E. This approach is favored in semiarid ecosystems in which there is a close coupling of ET and GPP and E makes up a considerable amount of monthly ET.

Several recently developed methods for partitioning eddy-covariance-measured ET are based on the optimality theory assumption that plants minimize water loss per unit of CO₂ gain (e.g., Hari et al., 2000; Katul et al., 2009; Medlyn et al., 2011; Schymanski et al., 2007). An outcome of this approach is that plant water use efficiency WUE, defined here as GPP∕T, scales with D^0.5 from which a relationship between GPP and T can be derived (Katul et al., 2009). Berkelhammer et al. (2016) noted that ET follows a linear relationship to GPP×D^0.5 and further assumed that the T ∕ ET ratio intermittently approaches 1. They then separated ET measurements from eddy covariance into GPP classes for which a minimum ET, min(ET)|_GPP, can be defined. T ∕ ET can then be calculated using

$$\begin{array}{}\text{(3)}& T/\mathrm{ET}={\displaystyle \frac{\mathrm{ET}}{min\left(\mathrm{ET}\right){\mathrm{|}}_{\mathrm{GPP}}}}.\end{array}$$

Applying this approach to different forests revealed considerable synoptic-scale variability in T ∕ ET that was dampened at seasonal timescales and compared well against isotopic approaches (Berkelhammer et al., 2016).

Zhou et al. (2016) built upon earlier work (Zhou et al., 2014) and assumed that an ecosystem has an actual underlying water use efficiency (uWUE_a, where WUE in this case is defined as GPP∕ET), which is maximal or reaches its potential underlying water use efficiency (uWUE_p) when T ∕ ET approaches unity. T ∕ ET can thus be calculated from the ratio of actual to potential uWUE using optimality assumptions for both:

$$\begin{array}{}\text{(4)}& {\mathrm{uWUE}}_{\mathrm{p}}={\displaystyle \frac{\mathrm{GPP}\sqrt{D}}{T}}\end{array}$$

and

$$\begin{array}{}\text{(5)}& {\mathrm{uWUE}}_{\mathrm{a}}={\displaystyle \frac{\mathrm{GPP}\sqrt{D}}{\mathrm{ET}}}.\end{array}$$

Again, assuming that T ∕ ET intermittently approaches 1 in sub-daily eddy covariance measurements, the uWUE_p can be estimated empirically using 95th quantile regression to find the upper boundary of the relationship between measured ET and GPP×D^0.5. uWUE_a can be calculated using eddy covariance observations, and T estimates using this approach compare well against independent sap flow measurements (Zhou et al., 2018) and expected responses to drought (Han et al., 2018). A semiempirical model based on the uWUE concept by Boese et al. (2017) included radiation and was able to outperform the Zhou et al. (2016) approach, on average, consistent with the notion that T is also driven by radiation (Eq. 1) (Pieruschka et al., 2010). It is important to note when applying WUE-based approaches that there are important discrepancies between WUE measurements at the leaf and canopy scales that still need to be resolved (Medlyn et al., 2017; Medrano et al., 2015) and also that GPP estimates from eddy covariance observations may have considerable uncertainty.

In a more sophisticated attempt to partition ET utilizing optimality theory, Perez-Priego et al. (2018) utilized a big-leaf canopy model in which parameters were optimized using half-hourly data in 5-day windows. Uniquely, the marginal carbon cost of water was factored into the cost function during parameter estimation, so the parameters for each 5-day window maximized the fit between modeled and observed GPP and also minimized water loss per carbon gain. T was then calculated using g_c from the model, and E was calculated as the residual (ET−T).

A modified (in this case binned) parameter optimization approach was used by Li et al. (2019) to estimate g_surf, which follows the model proposed by Lin et al. (2018):

$$\begin{array}{}\text{(6)}& g_{\mathrm{surf}}=g_{\mathrm{0}}+g_{\mathrm{1}}{\displaystyle \frac{\mathrm{GPP}}{D_{\mathrm{L}}^m}}.\end{array}$$

Here, g₀ (assumed to correspond to soil conductance), g₁ (assumed to correspond to vegetation conductance), and m are optimized parameters, D_L is the inferred leaf-level D, and g_surf is estimated by inverting Eq. (1) and is assumed to represent ecosystem conductance to water vapor flux. Rather than optimizing using a moving window over time, data were binned using independent soil moisture data associated with the eddy covariance site, with g₀, g₁, and m optimized in each bin to account for changes due to water limitations. Partitioning was then calculated as

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{T}{\mathrm{ET}}}={\displaystyle \frac{g_{\mathrm{1}}}{g_{\mathrm{surf}}}}\end{array}$$

and

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{E}{\mathrm{ET}}}={\displaystyle \frac{g_{\mathrm{0}}}{g_{\mathrm{surf}}}}.\end{array}$$

The Perez-Priego et al. (2018) and Li et al. (2019) methods both circumvent the assumption that T ∕ ET approaches unity at some periods by estimating ecosystem conductances directly. The transpiration estimation algorithm (TEA) from Nelson et al. (2018) utilizes a nonparametric model and thereby further limits assumptions made about how the ecosystem functions. However, TEA must make the assumption that T ∕ ET approaches 1, which it does by removing observations when the surface is likely to be wet. In a validation study that utilized model output as synthetic eddy covariance datasets in which E and T are known, TEA was able to predict T ∕ ET patterns in both space and time but showed a sensitivity to the minimum modeled E. Overall, TEA was able to predict temporal patterns of T across three different ecosystem models and provides an important basis for comparison because the model for T is agnostic to underlying ecosystem function.

3.2 Partitioning ET using high-frequency eddy covariance observations

Scanlon and Kustas (2010) (see also Scanlon and Sahu, 2008) developed a partitioning approach for E and T using high-frequency eddy covariance measurements based on the notion that atmospheric eddies transporting CO₂ and water vapor from stomatal processes (T and net primary production; NPP=GPP – aboveground respiration by the autotrophic canopy) and non-stomatal processes (E and soil respiration) independently follow flux–variance similarity as predicted by Monin–Obukhov similarity theory. In brief, there are two end-member scenarios for a parcel of air transported from a surface: one without stomata and one with stomata. An eddy transported away from a surface that is respiring CO₂ and evaporating water through pathways other than stomata will have deviations from the mean CO₂ mixing ratio (c^′) and water vapor mixing ratio (q^′) that are positively correlated. An eddy of air transported by a surface with stomata will have a negative relationship between c^′ and q^′ due to CO₂ uptake and T during daytime, whose ratio can also be described by a unique WUE at the leaf level. This leaf-level WUE is thereby used to establish a functional relationship between the variance of CO₂ due to stomatal uptake (${\mathit{\sigma }}_{\mathrm{cp}}^{\mathrm{2}}$) and the correlation between stomatal and non-stomatal CO₂ exchange processes (ρ_cp,cr). Subsequently, ET can be partitioned into its T and E components by matching the observed correlation of q^′ and c^′ (ρ_q,c) to the corresponding value of ρ_cp,cr (Scanlon and Sahu, 2008). The original approach applied wavelet filtering to remove large-scale atmospheric effects that impact the validity of underlying flux–variance relationships and was shown to realistically reproduce T ∕ ET relationships over the growing period of a corn (maize) crop (Scanlon and Kustas, 2012).

Subsequent work by Skaggs et al. (2018) noted that there is an algebraic solution to terms that had previously been solved using optimization (namely ${\mathit{\sigma }}_{\mathrm{cp}}^{\mathrm{2}}$ and ρ_cp,cr; Palatella et al., 2014) and created an open-source Python module, fluxpart, to calculate E and T using the flux–variance similarity approach. The first applications of the flux–variance similarity approach used a leaf-level WUE formulation following Campbell and Norman (1998); fluxpart allows leaf-level WUE to vary as a function of D or take a constant value. Leaf-level WUE varies throughout the canopy and in response to other environmental conditions. Using high-frequency measurements above the canopy rather than leaf-level observations to estimate it results in uncertainties (Perez-Priego et al., 2018). These uncertainties in leaf-level WUE can be addressed in part by using outgoing longwave radiative flux density observations to estimate canopy temperature (Klosterhalfen et al., 2019a, b). A careful comparison of flux–variance partitioning results against fluxes simulated by large eddy simulation revealed that it yields better results with a developed plant canopy with a clear separation of CO₂ and water vapor sources and sinks (Klosterhalfen et al., 2019b). It is also possible to separate E and T using conditional sampling of turbulent eddies (Thomas et al., 2008); the performance of the conditional sampling method is a function of canopy height and leaf area index, and the performance of the flux–variance similarity method is related to the ratio between sensor height and canopy height (Klosterhalfen et al., 2019a), suggesting that different methods may deliver better results in different ecosystems with differing measurement setups.

It should also be noted that flux–variance similarity can be used directly with (half-)hourly flux data if the wavelet filtering step is negligible (necessary variables of each time period are the CO₂ and water vapor flux, their respective variances (${\mathit{\sigma }}_c^{\mathrm{2}}$, ${\mathit{\sigma }}_q^{\mathrm{2}})$, ρ_q,c, and an estimated leaf-level WUE), but in practice high-frequency eddy covariance data are required because the necessary terms are rarely computed and saved. Of course, all eddy-covariance-based ET partitioning approaches need to (i) take decoupling between atmosphere, canopy, and subcanopy into account (e.g., Jocher et al., 2017); (ii) critique the energy balance closure of the observations (Leuning et al., 2012; Stoy et al., 2013; Wohlfahrt et al., 2009), especially in closed-path eddy covariance systems that are prone to water vapor attenuation in the inlet tube (Fratini et al., 2012; Mammarella et al., 2009); and (iii) acknowledge the uncertainty of eddy-covariance-based GPP estimates. An advantage of eddy-covariance-based approaches to partition E and T is that they can be complemented by other new approaches that measure or estimate E and T at temporal scales that align with the common half-hourly or hourly eddy covariance averaging period and spatial scales that align with the eddy covariance flux footprint.

3.3 Solar-induced fluorescence (SIF)

GPP and T are coupled through stomatal function, and studies of GPP have recently been revolutionized by space- and ground-based observations of solar-induced fluorescence (SIF) (Frankenberg et al., 2011; Gu et al., 2018; Köhler et al., 2018; Meroni et al., 2009), the process by which some of the incoming radiation that is absorbed by the leaf is reemitted by chlorophyll. SIF emission is related to the light reactions of photosynthesis, but GPP estimation also requires information on the dark reactions and stomatal conductance such that the remote sensing community is currently challenged by how to use SIF to estimate GPP. New studies also propose that SIF might be used to monitor T, possibly in combination with surface temperature measurements, acknowledging the close link between GPP and T due to their joint dependence on stomatal conductance and common meteorological and environmental drivers (Alemohammad et al., 2017; Damm et al., 2018; Lu et al., 2018; Pagán et al., 2019; Shan et al., 2019).

While SIF is related to the electron transport rate (Zhang et al., 2014), T primarily depends on stomatal conductance such that SIF and T are linked empirically but not mechanistically. This link is expected if GPP and T are tightly coupled. SIF has also been proposed to predict the ecosystem-scale WUE (i.e., GPP∕T) (Lu et al., 2018), a critical component of many of the E and T partitioning algorithms based on the eddy covariance ET measurements described above. Shan et al. (2019) showed that T can be empirically derived from SIF in forest and crop ecosystems, with explained total variance ranging from 0.57 to 0.83, and to a lesser extent in grasslands with explained variance between 0.13 and 0.22. The authors suggested that the decoupling between GPP and T during water stress hampered the use of SIF to predict T, particularly in grasslands, noting that T can occur without GPP under periods of plant stress (Bunce, 1988; De Kauwe et al., 2019). There is a strong empirical link between the ratio of T over potential evaporation and the ratio of SIF over PAR, and the relationship depends on the atmospheric demand for water, with larger transpiration for the same SIF when potential evaporation is higher (Alemohammad et al., 2017; Damm et al., 2018; Lu et al., 2018; Pagán et al., 2019; Shan et al., 2019). These ratios vary with assumptions regarding the potential evaporation calculation as well (Fisher et al., 2010). SIF can be measured at multiple spatial and temporal scales (Köhler et al., 2018), including the scale of the eddy covariance flux footprint (Gu et al., 2018), and this information can in turn be incorporated into remote-sensing-based approaches for estimating ET using remote sensing platforms (see Appendix A) following additional mechanistic studies of its relationship with T.

3.4 Carbonyl sulfide (COS) flux

Other approaches to estimate GPP and g_c use independent tracers such as carbonyl sulfide (COS). When plants open their stomata to take up CO₂ for photosynthesis, they also take up COS (Campbell et al., 2008), a trace gas present in the atmosphere at a global average mole fraction of ∼500 ppt (Montzka et al., 2007). The leaf-scale uptake of COS, F_COS (pmol m⁻² s⁻¹), can be calculated using

$$\begin{array}{}\text{(9)}& F_{\mathrm{COS}}=-C_{\mathrm{COS}}{\left({\displaystyle \frac{\mathrm{1}}{g_{\mathrm{b}}}}+{\displaystyle \frac{\mathrm{1}}{g_{\mathrm{s},\mathrm{COS}}}}+{\displaystyle \frac{\mathrm{1}}{g_{\mathrm{i}}}}\right)}^{-\mathrm{1}},\end{array}$$

where C_COS (pmol mol⁻¹) is mole fraction of COS and g_b, g_s,COS, and g_i represent the leaf-scale boundary layer, stomatal, and internal conductances (here mol m⁻² s⁻¹) to COS exchange (Sandoval-Soto et al., 2005; Wohlfahrt et al., 2012). The latter lumps together the mesophyll conductance and the biochemical “conductance” imposed by the reaction rate of carbonic anhydrase, the enzyme ultimately responsible for the destruction of COS (Wehr et al., 2017). Equation (9) also makes the common assumption that, because the carbonic anhydrase is highly efficient in catalyzing COS, the COS mole fraction at the diffusion end point is effectively zero (Protoschill-Krebs et al., 1996). Provided appropriate vertical integration over the canopy is made, Eq. (9) can be used to describe canopy-scale F_COS (Wehr et al., 2017).

Because COS and CO₂ share a similar diffusion pathway into leaves and because the leaf exchange of COS is generally unidirectional, COS has been suggested (Sandoval-Soto et al., 2005; Seibt et al., 2010; Wohlfahrt et al., 2012) and demonstrated (Wehr et al., 2017; Yang et al., 2018; Spielmann et al., 2019) to present an independent proxy for estimating GPP. Motivated by the common boundary layer and stomatal conductances, there has been recent interest in using measurements of the COS exchange to estimate the canopy stomatal conductance to water vapor and by extension T (Asaf et al., 2013; Wehr et al., 2017; Yang et al., 2018). Solving for g_s,COS from Eq. (9) requires measurements of F_COS (e.g., by means of eddy covariance; Gerdel et al., 2017) and C_COS, while g_b and g_i are typically estimated based on models.

With g_s (and by canopy scaling g_c) determined this way and an estimate of aerodynamic conductance (the canopy analog to the leaf boundary layer conductance; Eq. 1), T may be derived by multiplication with the canopy-integrated leaf-to-air water vapor gradient. The first and to date only study to attempt this was conducted by Wehr et al. (2017), who demonstrated excellent correspondence with g_c estimated from ET measurements in a temperate deciduous forest. While stomata dominated the limitation of the COS uptake during most of the day, co-limitation by the biochemical “conductance” imposed by carbonic anhydrase was observed around noon. This finding is consistent with leaf-level studies by Sun et al. (2018) and suggests that g_i in Eq. (9) may not generally be negligible, even though Yang et al. (2018) found the bulk surface conductance of COS (i.e., all conductance terms in Eq. 9 lumped together) to correspond well with the surface conductance for water vapor inferred from ET. As soils may both emit and take up COS, ecosystem-scale COS flux measurements need to account for any soil exchange, even though typically the soil contribution is small (Maseyk et al., 2014; Whelan et al., 2018). One notable exception for larger soil F_COS fluxes occurs in some agricultural systems (Whelan et al., 2016) due in part to the relationship of F_COS with soil nitrogen (Kaisermann et al., 2018). Clearly, further studies are required in order to establish whether the complexities of and uncertainties associated with inferring g_s from Eq. (9) and non-stomatal fluxes make COS observations a sensible independent alternative for estimating canopy T.

3.5 Advances in thermal imaging

Thermal remote sensing measures the radiometric surface temperature following the Stefan–Boltzmann law. ET can be estimated using thermal remote sensing by applying an ecosystem energy balance residual approach: $\mathit{\lambda }E=R_{\mathrm{n}}-G-H$ (Norman et al., 1995). Quantifying the available energy term (R_n−G) is difficult from space, and the radiometric surface temperature measured by infrared sensors is different from the aerodynamic surface temperature that gives rise to sensible heat flux (H) (Kustas and Norman, 1996). Despite these challenges, thermal remote sensing for ET has been widely used with multiple satellite platforms including Landsat, MODIS, Sentinel, and GOES (Anderson et al., 2012; Fisher et al., 2017; Semmens et al., 2016). One of NASA's newest missions is ECOSTRESS, mounted on the International Space Station, which produces thermally derived ET at 70 m resolution with diurnal sampling (Fisher et al., 2017).

Advances in thermal imaging (thermography) have made it possible to make radiometric surface temperature observations at increasingly fine spatial and temporal resolutions (Jones, 2004), on the order of millimeters or less, such that E and T can be measured individually from the surfaces from which they arise. Thermography has been used to estimate E_soil (Haghighi and Or, 2015; Nachshon et al., 2011; Shahraeeni and Or, 2010) and T from plant canopies (Jones, 1999; Jones et al., 2002), often in agricultural settings (Ishimwe et al., 2014; Vadivambal and Jayas, 2010). Researchers are increasingly using tower and UAV-mounted thermal cameras to measure the temperatures of different ecosystem components at high temporal and spatial resolution (Hoffmann et al., 2016; Pau et al., 2018), which could revolutionize the measurement of T from plant canopies (Aubrecht et al., 2016) or even individual leaves in a field setting (Page et al., 2018). Such measurements need to consider simultaneous E_i and T from wet leaf surfaces.

3.6 The challenges of measuring evaporation from canopy interception

E_i from wet canopies can return 15 %–30 % or more of incident precipitation back into the atmosphere annually (Crockford and Richardson, 2000), and models struggle to simulate it accurately (De Kauwe et al., 2013). Although interception has been studied for over a century, the underlying physical processes, atmospheric conditions, and canopy characteristics that affect it are poorly understood (van Dijk et al., 2015). Accurately estimating E_i from wet canopies is critical for the proper simulation of interception loss (Pereira et al., 2016). However, E_i predicted by the Penman–Monteith equation (Eq. 1) during rainfall is often a factor of 2 or more smaller than the E_i derived from canopy water budget measurements (Schellekens et al., 1999). A recent study using detailed meteorological measurements from a flux tower indicates that the underestimated E_i by the Penman–Monteith equation might be attributed to the failure in accounting for the downward sensible heat flux and heat release from canopy biomass, which can be major energy sources for wet-canopy E (Cisneros Vaca et al., 2018). Storm characteristics (e.g., amount, storm duration, and intensity) and canopy structural information (e.g., canopy openness, canopy storage capacity) are all important parameters for modeling E_i (van Dijk et al., 2015; Linhoss and Siegert, 2016; Wohlfahrt et al., 2006). To partition total ET into T, E_soil, and E_i, it is necessary to simulate the dynamics of canopy wetness before, during, and after each storm so that models can be applied to the dry and wet portions of the canopy, respectively (Liu et al., 1998), a process that can be implemented using a running canopy water balance model (Liu, 2001; Rutter et al., 1971; Wang et al., 2007). Understanding the sources of water is therefore useful for quantifying differences among T, E_soil, and E_i, and information from water isotopes can be helpful to do so.

3.7 Isotopic approaches

The hydrogen and oxygen atoms of water molecules exist in multiple isotopic forms, including ²H and ¹⁸O, which are stable in the environment and can be used to trace the movement of water through hydrologic pathways (Bowen and Good, 2015; Gat, 1996; Good et al., 2015; Kendall and McDonnell, 2012). Because heavier atoms preferentially remain in the more condensed form during phase change, evaporation enriches soils in ²H and ¹⁸O (Allison and Barnes, 1983), while root water uptake typically removes water from the soil without changing its isotope ratio (Flanagan and Ehleringer, 1991). This difference in the isotope ratio, $R=\left[{}^{\mathrm{2}}\mathrm{H}\right]/\left[{}^{\mathrm{1}}\mathrm{H}\right]$ or [¹⁸O] ∕ [¹⁶O], of E_soil compared with the isotope ratio of water moving through plants is the basis for the isotopic partitioning of ET. If ET consists of two components, E and T, with distinct isotopic composition – R_E for soil evaporation and R_T for plant transpiration – then the bulk flux, R_ET, can be incorporated into a simple mass balance of the rate isotope (i.e., $R_{\mathrm{ET}}\mathrm{ET}=R_{\mathrm{E}}E+R_{\mathrm{T}}T$), which can be rearranged as (Yakir and Sternberg, 2000)

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{T}{\mathrm{ET}}}={\displaystyle \frac{R_{\mathrm{ET}}-R_{\mathrm{E}}}{R_{\mathrm{T}}-R_{\mathrm{E}}}}.\end{array}$$

Thus, knowledge of the isotopic ratio of each flux component, R_E and R_T, as well as the total bulk flux isotope ratio, R_ET, is sufficient to estimate the fraction that passes through plants.

Techniques to measure R_ET have diversified since the widespread deployment of laser-based integrated cavity output spectroscopy (ICOS) systems, which currently monitor atmospheric stable isotope ratios, R_A, at a wide number of sites (Wei et al., 2019; Welp et al., 2012). Vertical profiles and high-frequency measurements of R_A are used to determine R_ET using multiple methods, all of which are associated with potentially large uncertainty (Griffis et al., 2005, 2010; Keeling, 1958). Propagation of uncertainties through Eq. (10) demonstrates that errors in R_ET, R_T, and R_E, as well as differences between R_E and R_T, strongly influence the final partitioning estimate (Good et al., 2014; Phillips and Gregg, 2001). The isotopic approach becomes uninformative as R_E approaches R_T. Furthermore, as E_i adds another source term to the isotope mass balance, Eq. (10) can be implemented over short periods only when the canopy is dry. If E_i is incorporated as a third source, its magnitude and isotope ratio must be specified, and these assumptions can strongly influence any final isotope-based partitioning estimates (Coenders-Gerrits et al., 2014; Schlesinger and Jasechko, 2014).

The value of R_E is derived from the soil water isotope ratio, R_S, as well as the temperature and humidity conditions under which evaporation happened (Craig and Gordon, 1965). The destructive extraction of water from soil cores can be used to estimate R_S, though recent studies have highlighted discrepancies between methodologies (Orlowski et al., 2016a, b). In situ monitoring of R_S obtained by pumping soil vapor through ICOS systems has been demonstrated (Gaj et al., 2016; Oerter et al., 2016; Volkmann and Weiler, 2014) and recently applied to ET partitioning to provide continuous updates on soil isotope ratios (Quade et al., 2019). Eddy covariance measurements of ²H and ¹⁸O are now possible (Braden-Behrens et al., 2019). However, identifying R_S remains challenging, and the bulk soil moisture composition (Mathieu and Bariac, 1996; Soderberg et al., 2013), depth (Braud et al., 2005), and soil physical composition (Oerter et al., 2014) at which evaporation occurs can alter the R_S to R_E relationship.

If water entering the plant is isotopically the same as transpired water, known as the isotopic steady-state assumption, then R_T=R_S. However, preferential uptake at the root–soil interface, differences between plant internal water pools in time, and mixing along the water pathways within plants will invalidate the steady-state assumption (Farquhar and Cernusak, 2005; Ogée et al., 2007). Finally, variability between and within plant species and plant–soil microclimates of an ecosystem will move the system away from the simple two-source model used in Eq. (10). Accurate knowledge of the isotope ratio within various water reservoirs of a landscape, including the planetary boundary layer (Noone et al., 2013), and how these translate into distinct water fluxes is required to advanced isotope-based partitioning approaches.

3.8 Statistical approaches

In addition to modeling g_surf as the sum of g_c, g_soil, and g_i, daily g_surf can also be well-approximated using emergent relationships between the atmospheric boundary layer and land surface fluxes, as demonstrated by the Evapotranspiration from Relative Humidity in Equilibrium (ETRHEQ) method (Rigden and Salvucci, 2015; Salvucci and Gentine, 2013). The ETRHEQ method is based on the hypothesis that the best-fit daily g_surf minimizes the vertical variance of relative humidity averaged over the day. Estimates of ET from this approach compare favorably to eddy covariance measurements (Gentine et al., 2016; Rigden and Salvucci, 2016), and the method can be applied at weather stations due to its primary dependence on meteorological observations. Rigden et al. (2018) recently developed a statistical approach to decompose estimates of g_surf from ETRHEQ into g_c and g_soil, allowing ET to be partitioned to T and E. The partitioning approach is based on the assumption that vegetation and soil respond independently to environmental variations and utilizes estimates of g_surf at ∼1600 US weather stations, meteorological observations, and satellite retrievals of soil moisture. Estimates of T from this statistical approach show strong agreement with SIF and realistic dry-down dynamics across the US (Rigden et al., 2018); however, the method lacks evaluation with E and T observations directly. Fortunately for many of the techniques discussed above, new large-scale methods for estimating E_soil based on theory have recently been developed and applied at large scales.

3.9 Novel approaches for estimating soil evaporation

E_soil is conventionally measured using lysimeters (Black et al., 1969), with some promising results from carefully designed chamber approaches that seek to minimize the impacts of the rapidly humidifying within-chamber atmosphere on evaporation (Raz-Yaseef et al., 2010; Yepez et al., 2005). E_soil has received extensive theoretical treatment (e.g., Brutsaert, 2014) that has resulted in models that align well with observations on ecosystem scales (e.g., Perez-Priego et al., 2018; Lehmann et al., 2018; Merlin et al., 2016, 2018). Lehmann et al. (2018) defined a new model for soil evaporative resistance that correctly describes the transition from stage-I evaporation (non-diffusion limited) to stage-II evaporation (diffusion limited). The model was able to correctly describe the soil moisture dependence of E_soil across different soil types. This approach was extended by Or and Lehmann (2019), who developed a conceptual model for soil evaporation called the surface evaporative capacitance (SEC) model for E_soil. Briefly, the transition between stage-I evaporation of a drying soil with capillary flow from deep moisture sources and stage-II evaporation characterized by water vapor diffusion is modeled using an evaporation characteristic length that differs by soil type (Lehmann et al., 2008, 2018). The SEC model accurately simulated E_soil datasets from different global regions, and adding global maps of precipitation and soil properties creates spatially distributed E_soil estimates to model global E_soil. The SEC model can be used in combination with other remotely sensed ET estimates (e.g., GLEAM; Appendix A) to partition ET.

4 Critiquing the assumptions of ET partitioning methods

4.1 Do ecosystems exhibit optimal responses to D?

Many WUE-based approaches for partitioning E and T (Sect. 3.1 and 3.2) hinge on the notion that g_c follows an optimal response to D. Recent data-driven studies have argued that g_c measured using eddy covariance is “slightly suboptimal”, averaging between D¹ and D^0.5 with a mean of D^0.55 rather than D^0.5 (Lin et al., 2018), or is “nearly optimal” and scales with GPP×D^0.55 (Zhou et al., 2015). Here, we test the assumption that plant canopies exhibit optimal responses to D by assuming that it serves as a constraint on WUE following an implication of optimality theory that 1 minus the ratio of leaf-internal CO₂ (c_i) to atmospheric CO₂ (c_a), $\mathrm{1}-c_{\mathrm{i}}/c_{\mathrm{a}}$, also scales with D^0.5 (see Eq. 18 in Katul et al., 2009). Using the definition of WUE as GPP ∕ T, expanding GPP and T using Fick's law, and excluding differences in mesophyll conductance,

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{GPP}}{T}}={\displaystyle \frac{g_{\mathrm{c}}\mathit{\epsilon }{c}_{\mathrm{a}}\left(\mathrm{1}-\frac{c_{\mathrm{i}}}{c_{\mathrm{a}}}\right)}{g_{\mathrm{c}}D}}.\end{array}$$

In this equation, g_c cancels and ε is the relative diffusivity of H₂O and CO₂ molecules. If ($\mathrm{1}-c_{\mathrm{i}}/c_{\mathrm{a}}$) scales with D^0.5, eddy-covariance-estimated WUE (i.e., GPP ∕ ET) should therefore scale with D^−0.5 if it can be assumed that measured ET approaches T. We tested this notion using micrometeorological and eddy covariance data from 240 sites that include ecosystem type and ecosystem energy balance measurements in the La Thuile FLUXNET database following Stoy et al. (2013). We assumed that E is a trivial component of ET when WUE values exceed 95 % of observations (Zhou et al., 2016) and use a boundary line analysis commonly used in studies of leaf and canopy conductance (Schäfer, 2011) to describe this 95 % threshold. We then took the mean of the upper 95 % of eddy covariance WUE observations in 0.3 kPa bins of D and fit an exponential model to these observations using nonlinear least squares (Fig. 4a) rather than fitting a linear model following log transformation for values that approach zero. Using this approach, we arrive at a mean (± standard deviation) exponential term of $-\mathrm{0.53}\pm \mathrm{0.17}$ from the 240 sites (Fig. 4b), which is not significantly different from −0.5 using a one-sample t test. Repeating this analysis with the FLUXNET2015 dataset reveals a mean exponential term of $-\mathrm{0.49}\pm \mathrm{0.15}$, which is likewise not different from −0.5.

Figure 4 (a) An example of a boundary line analysis to quantify the relationship between vapor pressure deficit (D) and ecosystem water use efficiency defined as gross primary productivity (GPP) divided by evapotranspiration (ET) for the case of a single ecosystem in the La Thuile FLUXNET database, in this case Vielsalm, Belgium (BE-Vie), using GPP∕ET values that represent the 95th percentile of 0.3 kPa D bins, shown as red dots. The model $\mathrm{GPP}/\mathrm{ET}=k{D}^m$, fit using nonlinear least squares, is shown as a red line. The value of m for BE-Vie is −0.53; the values $m=-\mathrm{0.5}$ following Medlyn et al. (2011) and $m=-\mathrm{1}$ following Leuning (1995) are shown for reference with the same fitted value of k. Individual half-hourly eddy covariance measurements are shown in light gray. (b) The distribution of the best-fit exponential parameter (m) for 240 sites in the La Thuile FLUXNET database that contained the full energy balance measurements and ecosystem type information used in Stoy et al. (2013).

Land surface models struggle to simulate this emergent property of ecosystems. Models for the ecosystems shown in Fig. 1 tend to dramatically overpredict the magnitude of the exponential term with a mean value of −2.9 (Table 2). The exponential term of the BEPS model was $-\mathrm{0.54}\pm \mathrm{0.06}$, similar to observations. Combined, these results suggest that an optimal canopy response to D may be a reasonable assumption despite the challenges of leaf-to-ecosystem scaling and despite the use of above-canopy rather than D_L here, but the considerable variability of the calculated exponential terms suggests that more research is necessary to understand conditions under which optimality is a reasonable assumption and when it is not. The discrepancy in calculated exponential terms between measurements and models further emphasizes the importance of improved carbon and water coupling in ecosystem models.

Table 2 The exponential term (m) of the model WUE = kD^m fit using nonlinear least squares to the 95th percentile of WUE values in 0.3 kPa bins of D (see Fig. 4). Sites: CA-Ca1 (Schwalm et al., 2007), CA-Obs (Griffis et al., 2003; Jarvis et al., 1997), US-Ho1 (Hollinger et al., 1999). Models: BEPS (Liu et al., 1999), CAN-IBIS (Williamson et al., 2008), CNCLASS (Arain et al., 2006), ECOSYS (Grant et al., 2005), ED2 (Medvigy et al., 2009), ISAM , ISOLSM (Riley et al., 2002), LOTEC (Hanson et al., 2004), ORCHIDEE (Krinner et al., 2005), SIB (Baker et al., 2008), SIBCASA (Schaefer et al., 2009), SSIB2 (Zhan et al., 2003), TECO (Weng and Luo, 2008). Data are available from Ricciuto et al. (2013).

4.2 Does T ∕ ET approach unity?

Also central to many E and T partitioning approaches is the notion that T ∕ ET intermittently approaches 1 (Berkelhammer et al., 2016; Nelson et al., 2018; Zhou et al., 2016; Wei et al., 2017), as suggested by modeling analyses and measurements (Wei et al., 2015, 2018). This assumption was critiqued by Perez-Priego et al. (2018), who demonstrated that T ∕ ET was rarely greater than 0.8 in a Mediterranean ecosystem, even during dry periods when surface soil moisture was less than 0.2 m³ m⁻³, and that E scaled with time^−0.5 following Brutsaert (2014) (see also Boese et al., 2019, and Li et al., 2019). These findings of a sustained evaporation component and nonzero E∕ET even during dry conditions were also supported by lysimeter measurements in a semiarid grassland (Moran et al., 2009) and partly confirmed by a recent study based on isotopes in shrubs and a steppe ecosystem (Wang et al., 2018). The maximum daily T ∕ ET found by Scanlon and Kustas (2012) in a maize agroecosystem was also about 0.8, but Rana et al. (2018) found daily values that intermittently exceeded 0.9 in wheat and fava bean fields, and multi-method comparisons suggest that T ∕ ET often approaches 0.85 (Rafi et al., 2019). Anderson et al. (2017a) found that T ∕ ET routinely exceeded 0.9 in sugarcane, with maximum daily values above 0.95, and Li et al. (2019) also found values greater than 0.9 for other crops. We can critique the notion that T ∕ ET approaches 1 by applying the flux–variance similarity partitioning approach to a wheat canopy from central Montana, USA, measured by Vick et al. (2016). Wheat has a characteristically high surface conductance (Bonan, 2008) and approaches an ideal transpiring surface during the main growth period (Bonan, 2008; Priestley and Taylor, 1972). The dryland wheat crops studied here draw water from depth such that surface soils are often dry (Vick et al., 2016), minimizing E_soil. Applying the flux–variance similarity method of Scanlon and Kustas (2010) to the wheat crop and allowing the algorithm to estimate water use efficiency suggests that T ∕ ET frequently exceeds 0.95 during daytime periods when the algorithm converges (Fig. 5a). Repeating this analysis for a winter wheat crop near Sun River, Montana, USA, using the flux–variance similarity algorithm of Skaggs et al. (2018) confirms this finding with an even higher proportion of T ∕ ET values (20 %) that exceed 0.95 when allowing the algorithm to estimate water use efficiency. T ∕ ET, however, exceeded 0.95 in less than 2 % of measurements using the approach of Perez-Priego et al. (2018) in a Mediterranean savanna ecosystem (Fig. 6). These observations suggest that the notion that T ∕ ET approaches 1 is a good assumption in some ecosystems, perhaps in ecosystems with high LAI, with implications for flux partitioning by the methods that rely on this assumption.

Figure 5 Distribution of the fraction of evapotranspiration arising during daytime hours from transpiration (T) and evaporation (E) estimated using the flux–variance similarity approach of Scanlon and Kustas (2010) from (a) a winter wheat field near Moore, MT, USA, described in Vick et al. (2016) using a version of the original algorithm. (b) A winter wheat field near Sun River, MT, using fluxpart (Skaggs et al., 2018).

Figure 6 The distribution of the T ∕ ET ratio for half-hourly observations from the partitioning approach of Perez-Priego et al. (2018) for the Majadas de Tietar (ES-Lma), Spain, research site.

5 Research imperatives

Few field experiments have sought to constrain ecosystem E and T estimates using multiple observations to quantify their response to environmental variability and to test the assumptions of partitioning approaches (Perez-Priego et al., 2017, 2018). Those that have note large discrepancies in T ∕ ET estimates from different techniques (Quade et al., 2019). Despite these challenges, a multi-measurement approach is necessary to understand different ecosystem water flux terms (Li et al., 2018), but most multi-method ecosystem-scale experiments using eddy covariance measurements seek to constrain the carbon cycle rather than the water cycle to which it is coupled (Hanson et al., 2004; Williams et al., 2009). Here, we outline the basics of an ecosystem-scale experiment designed to address uncertainties in E and T measurements (Fig. 7).

Figure 7 A schematic of an ecosystem experiment designed to measure transpiration and evaporation from soil and intercepted water using multiple complementary measurement approaches.

It would be best to introduce such an experiment in an ecosystem with a relatively simple species distribution and a clear separation of above- and below-canopy E and T sources to apply flux–variance approaches (Klosterhalfen et al., 2019b; Williams et al., 2004) before addressing more complex ecosystems with multiple canopy layers (Fu et al., 2018; Santos et al., 2016). Observations should occur on timescales commensurate with satellite remote sensing overpasses (see Appendix A); the half-hourly time step used in most eddy covariance observations is likely sufficient to approximate conditions captured by polar-orbiting satellites. For example, MODIS has a morning local overpass time for TERRA and afternoon overpass for AQUA; GOME2 makes SIF observations in the morning, and OCO-2 flies over in the early afternoon. There will be more opportunities to study diurnal patterns in E and T with the forthcoming and ongoing OCO-3, ECOSTRESS, and Geostationary Carbon Cycle Observatory (GeoCarb) missions. There are also under-explored opportunities to study ET partitioning using geostationary satellites like GOES (Bradley et al., 2010), which compromises a temporal resolution on the order of minutes with a spatial resolution from the distant geostationary orbit on the order of kilometers. These length scales may be well-captured by scintillometry, and ET partitioning approaches that employ scintillometry are largely lacking to date. A short measurement time step to align with satellite overpasses is possible for chambers, lysimeters, and sapflux measurements but not some isotopic approaches (Fig. 7). Critically, thermography, SIF, and COS flux can also be measured at these timescales. An ideal E and T partitioning experiment would make them both above and below plant canopies, in conjunction with below-canopy eddy covariance, to isolate E_soil (Fig. 7). For full water balance accounting, observations of drainage from the rooting zone using drainage lysimeters, soil moisture at multiple soil levels spanning the root zone, the flow of water down plant stems (stemflow), leaf wetness sensors, measurements of the amount of water held in plants themselves, and of course multiple precipitation gauges are required. Such a multi-measurement approach would also create an opportunity to compare the performance of emerging technologies like distributed temperature sensing from fiber-optic cables (Schilperoort et al., 2018), modeling cosmic ray neutron fields for soil water source estimation (Andreasen et al., 2016), and global navigation satellite system reflectometry (GNSS-R) for soil moisture estimation (Zribi et al., 2018). It remains difficult to assimilate E and T measurements into models using conventional data assimilation techniques because observations may contain substantial bias error yet still provide valuable information (Williams et al., 2009). Emerging approaches from machine learning in the Earth and environmental sciences may therefore be particularly useful for combining the best information from different measurement techniques into a mass- and energy-conserving model of the surface–atmosphere exchange of water (Reichstein et al., 2019). Regardless of the specifics of the multi-measurement approach for constraining E and T measurements, we advocate more investment into the study of ET, “green water”, given its central importance in provisioning resources to an increasingly resource-scarce planet (Schyns et al., 2019).

6 Conclusion

New measurement techniques and analytical approaches for partitioning E and T at the ecosystem scale provide critical opportunities to improve land surface models, remote sensing products, and ultimately our understanding of the global water cycle. Ecosystem-scale experiments that measure E and T using multiple approaches are needed to understand how E and T respond differently to climate variability and change across different global ecosystems and also to critique the assumptions made by ET partitioning approaches to improve their skill. By strengthening our focus on the water cycle in studies of coupled carbon and water fluxes, our understanding of the role of the land surface in the climate system can only improve.

Appendix A

For completeness we briefly describe common algorithms used by remote sensing platforms for estimating E, T, and ET, noting that additional approaches exist and are under development (El Masri et al., 2019). Many widely used algorithms including SEBAL (Bastiaanssen et al., 1998), METRIC (Allen et al., 2007; Su, 2002), and SEBS (Su, 2002) use an energy balance approach that does not explicitly seek to separate T from E but remain highly valuable for water resource management and hydrology.

A1 PT-JPL

The Priestley–Taylor Jet Propulsion Lab (PT-JPL) global ET remote sensing retrieval algorithm (Fisher et al., 2008) is based on the potential evapotranspiration (PET) formulation of Priestley and Taylor (1972), which replaces the adiabatic terms in Eq. (1) with a parameter, α_PT, that takes a value of 1.26 under ideal evaporating conditions:

$$\begin{array}{}\text{(A1)}& \mathrm{PET}={\mathit{\alpha }}_{\mathrm{PT}}{\displaystyle \frac{s\left(R_{\mathrm{n}}-G\right)}{s+\mathit{\gamma }}}.\end{array}$$

To reduce PET to ET, Fisher et al. (2008) introduced ecophysiological constraint functions (f functions, unitless multipliers between 0 and 1) following Jarvis (1976). These are based on D, relative humidity (RH), the normalized difference vegetation index (NDVI), and the soil-adjusted vegetation index (SAVI; Huete, 1988). PT-JPL calculates T, E_soil, and E_i explicitly using

$$\begin{array}{}\text{(A2)}& {\displaystyle }\mathrm{ET}=T+E_{\mathrm{soil}}+E_{\mathrm{i}},\text{(A3)}& {\displaystyle }T=\left(\mathrm{1}-f_{\mathrm{wet}}\right)f_{\mathrm{g}}{f}_{\mathrm{Temp}}{f}_{\mathrm{M}}{\mathit{\alpha }}_{\mathrm{PT}}{\displaystyle \frac{s\left(R_{\mathrm{nc}}-G\right)}{s+\mathit{\gamma }}},\text{(A4)}& {\displaystyle }{E}_{\mathrm{soil}}=\left(f_{\mathrm{wet}}+f_{\mathrm{SM}}\left(\mathrm{1}-f_{\mathrm{wet}}\right)\right){\mathit{\alpha }}_{\mathrm{PT}}{\displaystyle \frac{s\left(R_{\mathrm{ns}}-G\right)}{s+\mathit{\gamma }}},\text{(A5)}& {\displaystyle }{E}_{\mathrm{i}}=f_{\mathrm{wet}}{\mathit{\alpha }}_{\mathrm{PT}}{\displaystyle \frac{s\left(R_{\mathrm{nc}}-G\right)}{s+\mathit{\gamma }}},\end{array}$$

where f_wet is relative surface wetness (RH⁴), f_g is green canopy fraction, f_Temp is a plant temperature constraint, f_M is a plant moisture constraint, and f_SM is a soil moisture constraint. R_nc and R_ns are net radiation absorbed by the canopy and soil, respectively. PT-JPL has been tested against measured ET from hundreds of FLUXNET sites worldwide, with a monthly average r² of 0.90 across all sites and a slope ∕ bias of 1.07 using in situ data (Fisher et al., 2008, 2009). The PT-JPL model forms the core ET retrieval algorithm in the ECOSTRESS mission (Fisher et al., 2017; Hulley et al., 2017) onboard the International Space Station. New applications of the PT-JPL algorithm have included canopy indices derived from CubeSats (Aragon et al., 2018).

A2 PM-MOD16

The PM-MOD16 algorithm estimates ET at 8-day intervals at 1 km² pixels across the global terrestrial surface with MODIS observations following Mu et al. (2011) using Eq. (A2). The PM-MOD16 algorithm follows the Penman–Monteith model (Eq. 1) rather than the Priestley–Taylor model (Eq. A1) by modeling conductance terms (or resistance terms as the inverse of the conductance terms in Eq. 1) rather than including f functions as in the PT-JPL algorithm. It explains ca. 86 % of the variability in eddy-covariance-observed ET from 46 sites in North America (Mu et al., 2011).

A3 GLEAM

The Global Land Evaporation Amsterdam Model (GLEAM) (Miralles et al., 2011a, b) also uses a Priestley–Taylor approach (Eq. A1) to estimate ET. It employs the Gash (1979) analytical model for canopy rainfall interception and semiempirical stress functions that vary between 0 and 1 (similar to the f functions of the PT-JPL model) to reduce PET to T for canopies with different characteristics. For T, this stress function is calculated based on the content of water in vegetation and the root zone. The former is approximated based on microwave vegetation optical depth and the latter is calculated using a multilayer soil model driven by observations of precipitation and updated through the assimilation of microwave surface soil moisture. Validation studies against eddy covariance data at daily timescales show average correlations typically ranging 0.81–0.86 (Martens et al., 2017).

A4 DTD

The dual-temperature difference (DTD) model follows the notion that diurnal changes in air and radiometric surface temperatures are related to surface–atmosphere heat flux (Norman et al., 2000). It has since been applied to MODIS observations to estimate ET (Guzinski et al., 2013) and to partition E and T using a Priestly–Taylor scheme described in Song et al. (2018). DTD estimates of E and T compared well to estimates derived using the flux–variance similarity algorithm of Skaggs et al. (2018).

A5 ALEXI–DisALEXI

Atmosphere–Land Exchange Inverse (ALEXI) is a multi-scale surface energy balance modeling system building on the two-source energy balance (TSEB) land surface representation of Norman et al. (1995). The TSEB partitions the composite radiometric surface temperature, Temp_rad, into soil and canopy temperatures, Temp_s and Temp_c, based on the local vegetation cover fraction apparent at the thermal sensor view angle, f(θ):

$$\begin{array}{}\text{(A6)}& {\mathrm{Temp}}_{\mathrm{rad}}\cong f\left(\mathit{\theta }\right){\mathrm{Temp}}_{\mathrm{c}}+\left(\mathrm{1}-f\left(\mathit{\theta }\right)\right){\mathrm{Temp}}_{\mathrm{s}}.\end{array}$$

With information about Temp_rad, LAI, and radiative forcing, the TSEB evaluates the soil (subscript “s”) and canopy (“c”) energy budgets separately by computing system and component fluxes of R_n, H, and ET (i.e., ET from soil and canopy), as well as G. Because angular effects are incorporated in Eq. (A6), the TSEB can accommodate thermal data acquired at off-nadir viewing angles by geostationary satellites. The TSEB has a built-in mechanism for detecting thermal signatures of stress in the soil and canopy. An initial iteration assumes that T is occurring at a potential (non-moisture-limited) rate, while E_soil is computed as a residual to the system energy budget. If the vegetation is stressed and transpiring at significantly less than the potential rate, T will be overestimated and the residual E_soil will become negative. Condensation onto the soil is unlikely midday on clear days, and therefore E_soil<0 is considered a signature of system stress. Under such circumstances, T is iteratively down-regulated until E_soil=0, noting that this assumption has been challenged by recent observations in some ecosystems (Perez-Priego et al., 2018).

For regional-scale applications of the TSEB, air temperature boundary conditions are difficult to specify with adequate accuracy due to localized land–atmosphere feedback. To overcome this limitation, the TSEB has been coupled with an atmospheric boundary layer (ABL) model, thereby simulating land–atmosphere feedback internally. In ALEXI (Anderson, 1997; Anderson et al., 2007), the TSEB is applied at two times during the morning ABL growth phase (between sunrise and local noon) using thermal infrared observations from geostationary satellites. Energy closure over this interval is provided by a simple slab ABL model (McNaughton and Spriggs, 1986), which relates the rise in air temperature in the mixed layer to the time-integrated influx of H from the land surface. As a result, ALEXI uses only time-differential temperature signals, thereby minimizing flux errors due to absolute sensor calibration and atmospheric and emissivity corrections (Kustas et al., 2001). For local-scale applications on length scales similar to many flux footprints on the order of 100 m, the coarse-scale flux estimates can be spatially disaggregated using the DisALEXI technique (Norman et al., 2003). DisALEXI uses air temperature diagnosed by ALEXI at a nominal blending height along with high-resolution LAI and land surface temperature data from polar-orbiting satellites to estimate fluxes at finer scales.