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

Abstract. 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.



Introduction
Some 70 000 km 3 of water leaves terrestrial ecosystems and enters the atmosphere through evapotranspiration (ET) every year Oki and Kanae, 2006). Despite its importance, we are unsure whether global ET has been increasing over time (Brutsaert, 2013(Brutsaert, , 2017Brutsaert 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).
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 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.

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) 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 . Upscaling this relationship results in a global estimate of terrestrial annual T / ET of 0.57±0.07 Biogeosciences, 16, 3747-3775, 2019 www.biogeosciences.net/16/3747/2019/ 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 , SIBCASA (Schaefer et al., 2009), SSIB2 (Zhan et al., 2003), TECO (Weng and Luo, 2008). Data are available from Ricciuto et al. (2013). (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 . 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 . 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 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 . 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).

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  Or and Lehman (2019) 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(Katul et al., , 2012. These complex dynamics of water storage and transport impact the conductance of water between ecosystems and the atmosphere (Mencuccini et al., 2019;, 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 −1 ) responds to represent the mass and energy balance of water flux between the land surface and the atmosphere (Monteith, 1965;Penman, 1948): In the Penman-Monteith equation, λ is the latent heat of vaporization (J kg −1 ), ρ is the density of water (kg m −3 ), s is the slope of the saturation vapor pressure function (Pa K −1 ), R n is the surface net radiation (W m −2 ), G is the ground heat flux (W m −2 ), ρ a is dry air density (kg m −3 ), c p is the specific heat capacity of air (J kg −1 K −1 ), D is the vapor pressure deficit (Pa), γ is the psychrometric constant (Pa K −1 ), g a is the conductance of the atmosphere, and g surf is surface conductance to water vapor flux (both m s −1 ). 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 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 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  (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 . 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 , but D is also highly coupled to soil moisture , and both depend on ET itself through soil-vegetationatmosphere coupling. Increases in atmospheric CO 2 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 2 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 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).
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.

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.

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. 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;. 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 abovecanopy airspace to accurately apply the eddy covariance technique in forested ecosystems (Jocher et al., 2017(Jocher et al., , 2018Paul-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 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 ET = m GPP + E 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 eddycovariance-measured ET are based on the optimality theory Biogeosciences, 16, 3747-3775, 2019 www.biogeosciences.net/16/3747/2019/ assumption that plants minimize water loss per unit of CO 2 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).  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 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: and 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  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 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 bigleaf 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): Here, g 0 (assumed to correspond to soil conductance), g 1 (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 0 , g 1 , and m optimized in each bin to account for changes due to water limitations. Partitioning was then calculated as and 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.

Partitioning ET using high-frequency eddy covariance observations
Scanlon and Kustas (2010)  ing high-frequency eddy covariance measurements based on the notion that atmospheric eddies transporting CO 2 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 2 and evaporating water through pathways other than stomata will have deviations from the mean CO 2 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 2 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 2 due to stomatal uptake (σ 2 cp ) and the correlation between stomatal and non-stomatal CO 2 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 fluxvariance 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 σ 2 cp 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 leaflevel 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 2 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 2 and water vapor flux, their respective variances (σ 2 c , σ 2 q ), ρ 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 halfhourly or hourly eddy covariance averaging period and spatial scales that align with the eddy covariance flux footprint.

Solar-induced fluorescence (SIF)
GPP and T are coupled through stomatal function, and studies of GPP have recently been revolutionized by spaceand 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 Biogeosciences, 16, 3747-3775, 2019 www.biogeosciences.net/16/3747/2019/ 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 .

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 2 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 −2 s −1 ), can be calculated using where C COS (pmol mol −1 ) 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 −2 s −1 ) 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 2 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 esti-mating 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 .

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: λE = R 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 Interna-tional Space Station, which produces thermally derived ET at 70 m resolution with diurnal sampling .
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.

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 wetcanopy E (Cisneros . 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.

Isotopic approaches
The hydrogen and oxygen atoms of water molecules exist in multiple isotopic forms, including 2 H and 18 O, which are stable in the environment and can be used to trace the movement of water through hydrologic pathways 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 2 H and 18 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, , 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 ET ET = R E E + R T T ), which can be rearranged as (Yakir and Sternberg, 2000) T 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(Griffis et al., , 2010Keeling, 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 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, Biogeosciences, 16, 3747-3775, 2019 www.biogeosciences.net/16/3747/2019/ 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 2 H and 18 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 twosource 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.

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 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.

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., 2016Merlin et al., , 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(Lehmann et al., , 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

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 1 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 2 (c i ) to atmospheric CO 2 (c a ), 1 − c i /c 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, In this equation, g c cancels and ε is the relative diffusivity of H 2 O and CO 2 molecules. If (1 − c i /c 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 −0.53 ± 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 −0.49 ± 0.15, which is likewise not different from −0.5. 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 −0.54 ± 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.

Does T / ET approach unity?
Also central to many E and T partitioning approaches is the notion that T / ET intermittently approaches 1 Nelson et al., 2018;Zhou et al., 2016;Wei et al., 2017), as suggested by modeling analyses and measurements (Wei et al., , 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  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). less than 0.2 m 3 m −3 , and that E scaled with time −0.5 following Brutsaert (2014) (see also Boese et al., 2019, andLi 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 . 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 fluxvariance 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 dur-   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), OR-CHIDEE (Krinner et al., 2005), SIB , SIBCASA (Schaefer et al., 2009), SSIB2 (Zhan et al., 2003), TECO (Weng and Luo, 2008). Data are available from Ricciuto et al. (2013).

Model
CACa1 CAObs USHo1 ing 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.

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 assump-  tions 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 , 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). 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 underexplored 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 wellcaptured 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  (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 . Regardless of the specifics of the multimeasurement 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).

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), MET-RIC (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: 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 where f wet is relative surface wetness (RH 4 ), 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 2 of 0.90 across all sites and a slope / bias of 1.07 using in situ data (Fisher et al., 2008(Fisher et al., , 2009). The PT-JPL model forms the core ET retrieval algorithm in the ECOSTRESS mission 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 2 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-covarianceobserved 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 multiscale 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 Temp rad ∼ = f (θ) Temp c + (1 − f (θ )) Temp s .
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 (nonmoisture-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.
Author contributions. PCS, TSEM, JBF, TG, OPP, and THS performed the analyses used to create figures, and PCS performed the analyses used to create Table 2. PCS, JBF, PG, TG, SPG, AK, SL, DGM, AJR, and GW were responsible for writing subsections of the paper, and all authors contributed to writing the text.
Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. Gabriel Bromley, Adam Cook, Elizabeth Vick, and Skylar Williams assisted with data collection and analysis, Chuck Merja, Scott Powell, James Irvine, and Bruce Maxwell provided logistical support, and Elke Eichelmann provided valuable comments on earlier drafts of this paper. We also acknowledge financial support for the eddy covariance data harmonization provided by CarboEuropeIP, FAO-GTOS-TCO, iLEAPS, the Max Planck Institute for Biogeochemistry, the National Science Foundation, the University of Tuscia, Universite´Laval, Environment Canada, the US Department of Energy, database development and technical support from the Berkeley Water Center, the Lawrence Berkeley National Laboratory, Microsoft Research eScience, Oak Ridge National Laboratory, the University of California -Berkeley, and the University of Virginia. We would additionally like to thank the North American Carbon Program Site-Level Interim Synthesis team, the Modeling and Synthesis Thematic Data Center, the Oak Ridge National Laboratory Distributed Active Archive Center, and Larry Flanagan for collecting, organizing, and distributing the model output used in this analysis. The authors want to thank Silvana Schott for her support in producing Fig. 7. Review statement. This paper was edited by Martin De Kauwe and reviewed by Zhongwang Wei and three anonymous referees.