Evapotranspiration regulates leaf temperature and respiration in dryland vegetation

Evapotranspiration regulates energy flux partitioning at the leaf surface, which in turn regulates leaf temperature. However, the mechanistic relationship between evapotranspiration and leaf temperature remains poorly constrained. In this study, we present a novel mechanistic model to predict leaf temperature as a linearized function of the evaporative fraction. The model is validated using measurements from infrared radiometers mounted on two flux towers in Arizona, USA, which measure canopies of Prosopis velutina with contrasting water availability. Both the observations and model predictions reveal that leaf temperature equilibrates with air temperature when latent heat flux consumes all of the energy incident on the leaf surface. Leaf temperature exceeds air temperature when there is a net input of energy into the leaf tissue. The flux tower observations revealed that evaporative cooling reduced canopy leaf temperature by ca. 1–5 ◦ C, depending on water availability. Evaporative cooling also enhanced net carbon uptake by reducing leaf respiration by ca. 15% in the middle of the growing season. The regulation of leaf temperature by evapotranspiration and the resulting impact on net carbon uptake represents an important link between plant water and carbon cycles that has received little attention in literature. The model presented here provides a mechanistic framework to quantify leaf evaporative cooling and examine its impacts on plant physiological function.


Introduction
Leaves serve as a critical nexus between water, energy, and carbon fluxes in terrestrial ecosystems, and leaf temperature (  ) plays an important role in regulating the rates of mass and energy fluxes at the leaf surface (Still et al., 2021;Vinod et al., 2022).  directly influences several physical processes that drive mass and energy exchange, including leaf-to-air vapor pressure deficit (VPD; Grossiord et al., 2020), thermal conductance and emittance (Jones, 2014), net photosynthetic assimilation (Medlyn et al., 2002), and leaf respiration (  ;Heskel et al., 2016).High values of   can also cause thermal stress and damage to leaf biochemical systems, which may permanently inhibit leaf physiologic function (O'Sullivan et al., 2017).  is thus a critical variable that regulates several aspects of terrestrial ecosystem function, and it is important to constrain the drivers of   to better predict the sensitivity of terrestrial ecosystems to anthropogenic climate change.
Generally speaking,   is regulated by environmental conditions and energy fluxes at the leaf surface.Empirical observations have demonstrated that   is often close to air temperature (  ), but the mechanistic relationship between   and   remains poorly constrained.Some studies have argued that leaves exhibit limited homeothermy, whereby the slope of the relationship between   and   is less than 1 (Michaletz et al., 2015(Michaletz et al., , 2016;;Blonder and Michaletz, 2018;Cook et al., 2021).Other studies have argued that leaves exhibit megathermy, whereby the slope of the relationship between   and   is greater than 1 (Salisbury and Spomer, 1964;Pau et al., 2018;Still et al., 2019bStill et al., , 2022)).Observations where   ≅   (i.e., poikilothermy) have also been reported (Drake et al., 2020;Miller et al., 2021;Uni et al., https://doi.org/10.1016/j.agrformet.2023.109560Received 11 November 2022; Received in revised form 16 March 2023; Accepted 6 June 2023 2022).The terminology for leaf thermal regimes follows the convention described by Cavaleri (2020).In practice,   observations are often normalized by   (i.e.,   −   ) to control for environmental variability and analyze other drivers of   .
Surface energy flux partitioning between latent () and sensible () heat flux also plays an important role in regulating   .Surface energy balance must be preserved at the leaf scale, so leaf-level  consumes energy that would otherwise increase   .Surface energy flux partitioning can be quantified using the evaporative fraction (  ), which measures the proportion of available energy (  ) that is consumed by : Thus, there is a direct physical relationship between   and   , which results in evaporative cooling of the leaf surface.
Evaporative cooling has important functional implications for plant carbon cycling and leaf physiologic function, particularly in hot and dry ecosystems (Hultine et al., 2020;Uni et al., 2022).Photosynthetic assimilation of carbon is highly dependent on   at the leaf scale (Medlyn et al., 2002).Maintaining lower   also reduces   (Heskel et al., 2016;Mathias and Trugman, 2022) and can prevent thermal damage to leaves (O'Sullivan et al., 2017).Indeed, several recent studies have speculated that plants may decouple photosynthesis and transpiration during extreme heat waves to maintain high levels of , which keeps   below critical thresholds that would result in damage to the leaf tissue (Drake et al., 2018;Krich et al., 2022;cf. De Kauwe et al., 2019).Likewise, water availability for evaporative cooling may limit the distributions of some plant species in dryland ecosystems when they cannot maintain physiologic function at ambient temperatures (Hultine et al., 2020).Improving mechanistic models of   will enhance our understanding of the feedbacks between water, energy, and carbon fluxes at the leaf surface and improve our ability to predict shifts in ecosystem function under anthropogenic climate change.It will also improve our ability to map ecosystem water fluxes at broad spatial scales using thermal remote sensing data (Mallick et al., 2022).
Many models predict   or   −  by combining energy balance theory with the Penman-Monteith equation (e.g., Monteith and Unsworth, 2013).However, implementing these models requires empirical assumptions about stomatal conductance, which is difficult to constrain.Here, we present an alternate modeling framework that predicts   as a linearized function of   .Our mechanistic model requires fewer surface parameters than previous formulations, which improves our ability to isolate and examine the environmental variables that drive   .The simplified model also yields fundamental insights into the relationship between   and   under varying environmental conditions, and the resulting impacts on plant physiologic function.
The paper is organized as follows.First, we present the new model.Then, we validate the model predictions using   measurements from infrared radiometers mounted on two flux towers in Arizona, USA, which measure stands of Prosopis velutina with contrasting water availability.We also examine the environmental variables that are most important for predicting   in the observational data set.Finally, we force the model with flux tower measurements to estimate the change in   that is attributable to evaporative cooling of the leaf surface, which may reveal an important link between evaporative cooling and net carbon uptake.In doing so, we address the following research questions: 1. How sensitive is   to changes in surface energy flux partitioning between  and ? 2. Which environmental variables directly regulate   ?Which of those variables is most important for regulating   in dryland ecosystems? 3. How much is   reduced by evaporative cooling of the leaf surface?

Leaf temperature model
Steady-state surface energy balance can be modeled as the difference between   , , and : When modeling the energy balance of individual leaves, the   term is equivalent to leaf-level net radiation (  ), which is the sum of downwelling (↓) and upwelling (↑) shortwave ( ) and longwave ( ) radiation fluxes: The  ↑ term can be calculated as a function of   measured in K: where   is leaf emissivity (  = 0.98),  is the Stefan-Boltzmann constant (5.67 ×10 -8 W m -2 K -4 ), and  is a coefficient that indicates whether it is a one-sided ( = 1) or two-sided ( = 2) leaf model.The  term in Eq. ( 2) can be calculated as: where  is air density,   is the specific heat of air, and   is the aerodynamic resistance to .We combined Eqs. ( 1), (2), and ( 5) and rearranged to produce a novel linearized equation for   −   : The expanded version of Eq. ( 6) contains two additive terms: a radiative heating term that is proportional to   and an evaporative cooling term that is proportional to   .Importantly, Eq. ( 6) reveals that   −   is a linear function of   and that the slope and intercept are functions of   and   .Because the slope and intercept negate each other (i.e., slope = intercept × -1),   −   = 0 • C when   = 1, which reveals that   converges at   when  consumes all of the energy incident on the leaf surface.It follows that: Eq. ( 7) provides a framework to examine the competing roles of   , radiative heating, and evaporative cooling in regulating   .It also provides a framework to examine the resulting impacts on plant physiologic function.
It is worth noting that   appears on both sides of Eq. ( 7) because   regulates   and thus   (i.e., Eq. ( 4)).If   is not measured directly, Eq. ( 7) can be solved numerically, as is discussed below in Section 2.3.7.Alternatively, the   term can be approximated using isothermal net radiation ( , ) following Jones (2014): We used direct measurements of   throughout our analysis, except as noted in Section 2.3.7.

Leaf temperature sensitivity analysis
We modeled the sensitivity of   −   to environmental drivers by forcing Eq. ( 6) with simulated values of   (250, 500, and 750 Wm -2 ),   (1, 10, 20, 30, and 40 sm -1 ), and   (0-1).The  term was held constant at 1.006 kg m -3 , and the   term was held constant at 1010 J K -1 kg -1 , which are representative values for the study area.A one-sided model ( = 1) was used to facilitate intercomparison with subsequent analyses.

Flux tower observations
We compared the modeled sensitivities from Eq. ( 6) to measurements from two eddy covariance flux towers in Arizona, USA, where   was measured by infrared radiometers mounted on the towers.The infrared radiometers measured the average temperature of many leaves on the outside of the P. velutina canopies, so we use the term   to describe the radiometer measurements of ''canopy-scale leaf temperature'' following (Still et al., 2021).However, we generally assume that   ≅   , and we use   and   interchangeably.The flux tower data set contains 17 site-years of growing season measurements under varying environmental conditions.We also used the flux tower data to analyze the environmental drivers of   and to quantify the impact of   on   .
The   measurements by the infrared radiometers represent leaves on the outside of a single P. velutina canopy, while the eddy covariance measurements represent the average fluxes within the fetch of the sensors (ca.50-200 m).We acknowledge the scale mismatch between the canopy-scale   measurements and the fetch-scale flux measurements, but we contend that novel insights can still be gleaned from the measurements using a ''big leaf'' assumption, whereby the entire fetch is assumed to behave like a single leaf in order to link leafscale theoretical models with canopy-scale measurements (e.g., Sellers et al., 1992;Amthor, 1994).In this context, Eq. ( 7) can be used to predict the average surface temperature of an entire canopy or stand.Canopy-scale processes are arguably more important than leaf-scale processes for understanding terrestrial ecosystem function, but they are also more difficult to constrain (Bonan, 2016).The temperature and fluxes of individual leaves can be measured using in situ sensors, but they may not be representative of the canopy as a whole (Miller et al., 2021;Vinod et al., 2022).We believe that the big leaf assumption is a reasonable approach to glean insights into the theoretical drivers of canopy-scale processes.To connect the leaf-scale theoretical model with the canopy-scale measurements, we included a ground heat flux () term when calculating   , such that: Eq. ( 9) helps control for the loss of available energy through the bottom of the canopy.We also used a one-sided model ( = 1) for all analyses.

Study sites
We analyzed data from two flux towers located in stands of P. velutina in southeastern Arizona, USA.Southeastern Arizona has a semi-arid climate with monsoonal precipitation that is delivered in brief, spatially restricted storms that dominate total annual rainfall and runoff (Thomas and Pool, 2006;Singer and Michaelides, 2017).The summer growing season encompasses both the driest and wettest parts of the year.The first part of the growing season is very dry, but ecosystems receive intense precipitation after the onset of the monsoon around early July.The monsoonal precipitation and accompanying humidity typically decrease after August, but generally remain above pre-monsoon levels through the end of the growing season (Higgins et al., 1997).The two stands that we analyzed have contrasting physiographic positions resulting in differences in plant water availability, particularly during the dry months before the onset of the monsoon.The differences in water availability create ideal conditions for a natural experiment to quantify the sensitivity of   to  and environmental conditions, while holding regional climatic variables relatively constant.
The first flux tower is in a riparian woodland approximately 16 km northeast of Sierra Vista, Arizona (31.6637 • N, 110.1777 • W).The riparian woodland is located on an old alluvial terrace above the San Pedro River, where the depth to groundwater is approximately 10 m (Sabathier et al., 2021).The flux tower is located ca. 225 m from the river channel, and the alluvial terrace is ca. 10 m above the river channel, so we assume that evaporation from the river channel does not affect the E measurements.The mean summer air temperature is 25 • C and the mean annual precipitation is 319 mm (PRISM Climate Group, 2021;Huntington et al., 2017).The woodland is dominated by a canopy of P. velutina (canopy cover ∼70%) with a mean height of 7 m and a maximum height of 10 m.Leaf emergence for the deciduous P. velutina trees typically occurs in April, and plant hydraulic function increases in late May.The understory is dominated by the perennial grass Sporobolus wrightii but annual forbs and herbs are common during the summer monsoon season (Scott et al., 2004).Rooting depths of P. velutina can exceed 10 m (Stromberg, 2013), and the pre-monsoon fluxes reveal that overstory vegetation accesses groundwater.Because groundwater provides a stable source of water that is somewhat decoupled from the local precipitation regime, evapotranspiration ( ) consistently exceeds precipitation on an annual basis (Missik et al., 2021;Scott et al., 2021).Groundwater is an important water source for maintaining vegetation structure and function in many dryland riparian plant communities (Kibler et al., 2021).The understory vegetation has a maximum rooting depth of 2-3 m, so it does not have access to groundwater and is dependent on water inputs from local precipitation (Scott et al., 2004).
The second flux tower is in an upland savanna at the Santa Rita Experimental Range, approximately 45 km south of Tucson, Arizona (31.8214 • N, 110.8661 • W).The site is a semi-desert grassland that has been encroached by P. velutina.The mean summer air temperature is 26 • C and the mean annual precipitation is 368 mm (PRISM Climate Group, 2021;Huntington et al., 2017).Scott et al. (2009) reported that the P. velutina canopy ranges in height from 0.25 to 6 m (mean height 2.5 m) and covers ∼35% of the ground area.Leaf emergence for P. velutina typically occurs in April (Seyednasrollah et al., 2019).The P. velutina plants at the upland savanna likely had lower leaf area index and smaller average leaf size than those at the riparian woodland (Stromberg et al., 1993).Perennial grasses, forbs, and subshrubs cover ∼22% of the ground area (Scott et al., 2009).Depth to groundwater exceeds 100 m, so the overstory and understory vegetation do not have access to groundwater and are dependent on water inputs from local precipitation.
Both flux towers contain an array of eddy covariance, meteorological, and soil sensors, along with infrared radiometers (IRT-P, Apogee Instruments, Logan, UT) pointed 45 • off-nadir at the P. velutina canopies.This study primarily relied on flux measurements of , ,   , , SW↓, SW↑, and LW↓; meteorological measurements of   , wind speed (), and relative humidity (RH); measurements of soil temperature (  ) and soil water content (SWC); and   measurements from the infrared radiometers.In the riparian woodland,  and  were measured at 14 m. was quantified for the surface using soil heat flux plate measurements at 5 cm depth along with the change in heat storage from 0-5 cm depth.Canopy-level   , , and RH were measured at 8 m.   was measured at 5 cm depth and SWC was measured at 22.5 cm depth.The infrared radiometer was mounted at 10 m.A four-component net radiometer measured individual SW↓, SW↑, LW↑, and LW↓ fluxes from 2001-2003, but it was replaced by a two-channel SW and LW net radiometer from 2004-2006.The four-channel radiometer was mounted at 14 m, and the sensors for the two-channel radiometer were mounted at 10 m and 14 m.In the upland savanna,  and  were measured at 7.8 m; SW↓, SW↑, LW↑, and LW↓ were measured at 7.1 m; and  was quantified for the surface.Canopy-level   and RH were measured at 2 m.Canopylevel  was measured at 3.5 m.   was measured at 5 cm depth and SWC was measured at 20 cm depth.The infrared radiometer was mounted at 7 m.Atmospheric transmittance and emittance were assumed to have a negligible impact on the radiometer measurements over the short distances (ca. 5 m) between the radiometers and the canopies (Aubrecht et al., 2016).We also assumed that the differences in   between the meteorological sensors and the canopies were negligible, given that the   sensors were at approximately the same heights as the measured leaves.All flux tower data were acquired from Ameriflux (sites US-CMW and US-SRM, respectively).See Scott (2021a,b) and Scott et al. (2004Scott et al. ( , 2009) ) for additional details about the collection and processing of the flux tower data.
Several criteria were used to filter the half-hourly flux tower measurements: • growing season observations between May and September • daytime observations between 8:00 and 16:00 local time • removed days with any measured precipitation and the day after any measured precipitation • observations with friction velocity ( * ) > 0.2 We also removed years that did not have complete records of growing season measurements and years where there were apparent shifts in the infrared radiometer view angle due to a loose mounting bracket, as evidenced by sudden changes in the relationship between   and   at 5 cm depth.Based on these criteria, two site-years of data were removed for the riparian woodland (2007 and 2008), and four site-years of data were removed for the upland savanna (2014, 2015, 2020, and 2021).The resulting data set contained six site-years of data for the riparian woodland (2001)(2002)(2003)(2004)(2005)(2006) and eleven site-years of data for the upland savanna (2007-2019, excluding 2014 and 2015).
The   analysis relied on individual measurements of SW↓, SW↑, and LW↓ fluxes, which were only available from 2001-2003 at the riparian woodland.They were available for all years at the upland savanna.The flux tower measurements were used to force Eq. ( 7) and analyze the sensitivity of   to  and environmental conditions.All of the other terms in Eq. ( 7) can be directly derived from the flux tower measurements, except for   .

Resistance to sensible heat flux
Following Young et al. (2021), the canopy-scale   is the sum of the resistance to momentum transfer (  ) and the excess resistance ( ℎ ): The   term can be estimated as a function of  and the friction velocity ( * ): The  ℎ term is a function of the roughness lengths for momentum ( 0 ) and heat ( 0ℎ ) as well as stability functions for momentum (  ) and heat ( ℎ ) exchange: where  is the Von Kármá n constant ( = 0.41).Eq. ( 8) can be simplified to ignore the stability functions, which have a negligible impact on the predicted values of  ℎ at a canopy scale (Young et al., 2021): The  0 and  0ℎ terms are often represented by the parameter  −1 such that: At an ecosystem scale, the parameter  −1 varies as a function of land cover, leaf area, vegetation structure, and environmental conditions (Yang and Friedl, 2003).Various empirical formulations for  −1 have been developed.We estimated  −1 as an empirical function of  * following (Thom, 1972), which yielded the most parsimonious predictions of   out of 12 formulas described by Verhoef et al. (1997) and Hong et al. (2012).The comparison of the formulas is described in the Supplementary Materials.

Model validation
The flux tower measurements were used to validate the model described in Eq. ( 7).We compared the   measurements from the infrared radiometers to   predictions that were generated by forcing Eq. ( 7) with concurrent flux tower measurements.The MAE, R 2 , slope, and bias were used to quantify the model performance at each site.

Energy balance closure
The model presented in Eq. ( 7) assumes energy balance closure.However, energy balance closure is rarely achieved in eddy covariance measurements due to systematic sensor errors, differences in the spatial footprints of individual sensors, advective fluxes, and a variety of other factors (Stoy et al., 2013;Mauder et al., 2020).The energy balance closure ratio () can be calculated as: We calculated  for the half-hourly flux measurements using Eq. ( 17).
We quantified the sensitivity of the   predictions to  by forcing closure in the flux measurements and then comparing the   predictions from the forced and unforced values.While forcing energy balance closure is often not recommended for eddy covariance analyses (e.g., Scott, 2010), comparing the different   predictions enabled us to quantify the model error that might be attributable to the lack of energy balance closure.Energy balance closure was forced by assuming that the  and  were measured correctly and adjusting the value of   .Energy balance closure can also be forced by assuming that   was measured correctly and adjusting the values of  and  (Twine et al., 2000;Knauer et al., 2018).However, the  term is not explicitly represented in Eq. ( 7).Energy balance closure was forced by setting   equal to the sum of the turbulent fluxes: where the subscript  denotes that the value was adjusted to force energy balance closure.Eq. ( 7) was forced with  , to generate a new set of   predictions.All other model forcings remained unchanged.We compared the two sets of   predictions to estimate the model error that might be attributable to the lack of energy balance closure.

T/ET partitioning
The model presented in Eq. ( 7) also assumes that all  is attributable to leaf transpiration ( ).However, eddy covariance measurements are collected at a stand scale, and soil evaporation () may also contribute to the  signal.The ratio of  ∕ can be used to quantify the extent to which the  signal is attributable to  .Several methods have been proposed to partition  and  in eddy covariance measurements (Stoy et al., 2019).We reanalyzed data from Scott et al. (2021) and Nelson et al. (2020a,b), who partitioned data for the riparian woodland and upland savanna, respectively, using the method proposed by Nelson et al. (2018).We used their daily estimates of  and  to calculate  ∕ for the two sites.The  ∕ analysis included all days where there was at least one half-hourly measurement in the filtered eddy covariance data set and an estimate of  ∕ from the published data sets.The resulting data set covered years 2005-2006 for the riparian woodland and 2007-2013 for the upland savanna.

Analysis of flux tower measurements
We conducted several analyses to identify the mechanistic basis for the model behavior using the flux tower measurements.We compared the distributions of the   −   and   measurements and calculated the seasonal and diurnal climatology of   −   at each site.We also produced seasonal and diurnal climatologies for the individual drivers of   , including   ,   ,   , and   .Spearman rank correlation was used to quantify the sensitivity of   and   −  to the individual drivers.The data were grouped by month to assess seasonal changes in the variables that drive   and   −   .Spearman rank correlation was also used to quantify the sensitivity of   ,   , and   to environmental variables measured by the flux towers, which may have an indirect effect on   and   −   .The environmental variables include SW↓, VPD, , and soil water content (SWC).VPD was calculated using the flux tower measurements of   and RH following Allen et al. (1998).

Leaf respiration model
We also analyzed the sensitivity of daytime leaf respiration (  ) to changes in   caused by  variability.Leaf respiration is a complex biochemical process that varies as a function of leaf mass per area, leaf nitrogen and phosphorus concentrations, photosynthetic carboxylation capacity,   , and other variables (Atkin et al., 2015).Leaf respiration is also inhibited by sunlight during the daytime (Kok, 1948;Heskel et al., 2013).In practice,   is often estimated as an empirical function of   (Mathias and Trugman, 2022).We estimated leaf dark respiration ( , ) following Heskel et al. (2016): where   is the reference temperature and  , is measured in units of μmol CO 2 m -2 s -1 .The  , (  ) parameter was set to 1.
Eq. ( 20) helps control for the light inhibition of   , which is not represented in the  , estimates (Way et al., 2015).We estimated  ,ℎ using two temperature forcings:   and modeled canopy temperature with no evaporative cooling ( , ).We modeled   by forcing Eq. ( 7) with flux tower measurements.The   term was calculated using Eqs.( 10)-( 16).We used modeled rather than measured   values to control for any effect of changing sensor calibrations over the multi-annual time series.We modeled  , by setting  to 0 in Eq. ( 7).We also calculated   as a function of its component fluxes (Eqs.(3)-( 4)) to account for the effect of   on  ↑.When all terms in Eq. ( 7) are directly measured, the feedback between   and  ↑ is implicitly encoded in the flux measurements.However, when combining measured and forced flux values (i.e., by setting  to 0), the feedback must be explicitly specified in the analytical formulation.The  ↓ term must also be multiplied by   .The equation for  , can be written as: Eq. ( 21) estimates the temperature of a non-transpiring canopy.The terms of the equation were forced with flux tower measurements.
To make Eq.( 21) analytically tractable, we rewrote the equation in the form of a quartic function and solved for  , using a numerical solver (NumPy v1.23.3;Harris et al., 2020).Physically unreasonable values where modeled  , was less than modeled   were likely due to the lack of energy balance closure and were removed from the analysis.The difference between the estimates of  ,ℎ using   and  , revealed the marginal change in  ,ℎ that is attributable to   variability caused by  (i.e.,  ,ℎ ∕  ()).This framework enabled us to estimate the decrease in daytime   caused by evaporative cooling of the leaf surface (  ).

Leaf temperature model
The mechanistic model of   presented in Eq. ( 7) reveals that there are three drivers of   : (1)   , (2) a radiative heating term that is proportional to   , and (3) an evaporative cooling term that is proportional to   .The model predicts that   converges to   when  consumes all of the energy incident on the leaf surface, regardless of environmental conditions (Fig. 1).When   < 1,   −   also varies as a function of   and   .Importantly, the model predicts that   −   ≥ 0 • C under all conditions, although there are environmental conditions when   −   approaches 0 • C even though   < 1.Specifically, the model predicts that   −   = 0 • C when   = 0 Wm -2 or when   = 0 sm -1 , regardless of the value of   .The first condition often happens around dawn and dusk, while the second condition is unrealistic in real-world settings (Young et al., 2021).

Model validation
The model of   presented in Eq. ( 7) was validated by forcing Eq. ( 7) with flux tower measurements and comparing the   predictions to concurrent   measurements from infrared radiometers mounted on the flux towers.The   measurements represented the average temperature of many leaves on the outside of a P. velutina canopy, and we generally assumed that   ≅   .The model yielded strong fits at both study sites (Fig. 2).The predictions for the riparian woodland exhibited a stronger fit (MAE = 2.67 • C).The predictions for the upland savanna exhibited a slightly weaker fit (MAE = 3.42 • C), but the range of observed   values was also larger.The model tended to slightly overestimate   at both sites (mean bias = 2.53 • C and 1.55 • C, respectively).

Energy balance closure
We also assessed the impact of energy balance closure on the   predictions.The median energy balance closure ratio () in the riparian woodland was 0.86 with an interquartile range of [0.75, 0.98].The median  in the upland savanna was 0.83 with an interquartile range of [0.75, 0.92] (Supplementary Figure 4).Forcing energy balance closure by adjusting the   value reduced the   predictions by 1.06 • C in the riparian woodland and 1.64 • C in the upland savanna.The MAE between the two sets of   predictions was 1.38 • C in the riparian woodland and 1.8 • C in the upland savanna.The   predictions based on  , were more parsimonious than the predictions based on the unforced values when compared to the infrared radiometer measurements at both sites.The R 2 values were 0.89 and 0.7 for the riparian woodland and upland savanna, respectively.The analysis suggests that the overestimation of   seen in Fig. 2 is due, in part, to the lack of energy balance closure in the flux data.This is supported by comparing the model prediction error to the  values.In the riparian woodland, the average model prediction error was near 0 • C when  ≈ 1.When  decreased below 1, the model prediction error increased monotonically.In the upland savanna, the model prediction error also increased as  decreased below 1, although the upland savanna exhibited a less clear trend (Supplementary Figure 5).

T/ET partitioning
The analysis of  ∕ values indicated that the  signal was dominated by  and not  at both sites (Fig. 3).In the riparian woodland, the median daily  ∕ value ranged from 0.87 to 0.92 for each month.In the upland savanna, the median daily  ∕ value ranged from 0.72 to 0.8 for each month.The upland savanna exhibited more variability in  ∕ values, suggesting that  may have contributed more error to the model predictions at that site.The  ∕ values were relatively consistent throughout the growing season and did not exhibit any apparent seasonal trend at either site.6) calculated using different values of available energy (  ), evaporative fraction (  ), and resistance to sensible heat flux (  ).Air density () was held constant at 1.006 kg m -3 .The specific heat of air (  ) was held constant at 1010 J K -1 kg -1 .Fig. 2. Model predicted   compared to   measurements from the infrared radiometers for the riparian woodland (a) and upland savanna (b).The predicted   values were calculated by forcing Eq. ( 7) with flux tower measurements.The blue lines are the 1:1 line.The mean absolute error (MAE) is also indicated.White areas in the plots indicate that there were 0 observations in that portion of the feature space.The color scale saturates when there are more than 100 half-hourly observations in a given portion of the feature space.

Flux tower observations
We also analyzed the flux tower observations to characterize the mechanistic basis for the model behavior.In the riparian woodland,   −  remained close to 0 • C for the entire time series, although   −  varied both seasonally and diurnally.The mean   −   value across the entire data set was 1.57 • C (Fig. 4).  −  values near 0 • C indicate that there was a substantial degree of evaporative cooling in the riparian woodland.Otherwise,   would substantially exceed   because of energy inputs from solar radiation.The largest   −   values tended to occur in May when the trees were leafing out and plant hydraulic function was still increasing.The values decreased substantially starting in June (Fig. 5).In May, the average peak value was 4.27 • C, and by June the average peak value decreased to 2.44 • C, with lower peaks occurring in subsequent months.The daily maximum values tended to occur around 11:30 local time.
The upland savanna exhibited similar trends.The mean   −  value across the entire data set was 6.46 • C (Fig. 4).The largest values tended to occur in the dry months of May and June and decreased substantially starting in July when the summer rainy season began (Fig. 5).In June, the average peak value was 10.59 • C. By August, the average peak value decreased to 5.72 • C. Unlike the riparian woodland, the daily maximum values tended to occur around 13:00 local time.The sensitivity of   −   to   from the flux tower measurements was consistent with the sensitivity predicted by Eq. ( 6).At both sites,   −   converged to 0 • C as   approached 1 (Fig. 4), which is consistent with the model predictions shown in Fig. 1.When   < 1, the range and distribution of measured   −   values was also similar to the modeled values.The   −   values ranged from −1.24 • C (1 st percentile) to 7.61 • C (99th percentile) in the riparian woodland and −0.25 • C (1 st percentile) to 16.44 • C (99th percentile) in the upland savanna.The maximum values at both sites occurred when   approached 0, consistent with the model predictions.The   values for the upland savanna (mean   = 0.15) were on average lower than the   values for the riparian woodland (mean   = 0.46), which provides a mechanistic explanation for why   −   was generally greater at the upland savanna than the riparian woodland.
The flux tower measurements exhibited a non-linear relationship between   and   −  , especially at the riparian woodland.The model in Eq. ( 6) predicts a linear relationship between   and   −  when all other variables are held constant.The apparent non-linear relationship between   and   −   is likely due to covariance between   ,   , and   on seasonal and diurnal time scales.The model validation accounted for the changing values of   ,   , and   , and it demonstrated strong model performance at both sites (Fig. 2).

Drivers of leaf temperature
Seasonal climatologies of   ,   , and   revealed seasonal changes in the environmental variables that drive   −   .The   exhibited the most pronounced seasonal trends and generally tracked the onset of the monsoon.In the riparian woodland, median   increased from 0.2 in May to 0.57 in August.In the upland savanna, median   remained low in May and June (0.08 and 0.07, respectively) and increased to 0.32 in August.At both sites,   exhibited a much less pronounced seasonal trend.Monthly median   values ranged from 12.1 to 14.6 sm -1 in the riparian woodland and 20.4 to 24.6 sm -1 in the upland savanna.Likewise, median monthly   (measured at 12:00 local time) exhibited little seasonal trend and ranged from 604 to 652 Wm -2 in the riparian woodland and 501 to 528 -2 in the upland savanna.The greater values of   in the riparian woodland are likely due to the greater canopy cover with lower albedo as well as smaller  flux.The median were 9.3% and 15.2% in the riparian and upland savanna, respectively.The median  fluxes were 59 Wm -2 and 119 Wm -2 , respectively.The   term can also be calculated directly from temperature measurements by inverting Eq. ( 5) (Verhoef et al., 1997).The inversion method yielded a different seasonal trend, indicating that   decreased throughout the growing season.However, the values of   were generally similar using both methods (Supplementary Figure 7).
Spearman rank correlation was used to quantify the sensitivity of observed   to the individual variables that drive   , including   ,   ,   , and   .As expected,   was highly correlated with   in all months at both sites (  ≥ 0.75; Supplementary Figure 8).We controlled for   by repeating the analysis with   −   values.There were coherent seasonal trends in the correlations between   −   and   ,   , and   at both sites (Fig. 7).In the riparian woodland,   −   was highly correlated with   early in the growing season (  = −0.79 in May), but the sensitivity to   decreased as monsoonal moisture accumulated in the ecosystem (  = −0.14 in September).The sensitivity to   peaked in the middle of the summer (  = 0.71 in July) and was lower at the beginning and end of the growing season.The sensitivity to modeled   was negligible in all months (  ≤ 0.1).In the upland savanna,   −   was more sensitive to   in May and June (  = 0.56 and 0.58, respectively) and more sensitive to   after the onset of the monsoon in July.The sensitivity to   peaked in July (  = −0.68)and decreased at the end of the growing season.The sensitivity to modeled   was weak in all months (−0.15 ≤   ≤ 0.09).
Spearman rank correlation was also used to quantify the sensitivity of observed   ,   , and   to environmental variables measured by the flux towers, including SW↓, VPD, , and SWC.The   term was negatively correlated with SW↓ and VPD at both sites (Table 1).The   term was also negatively correlated with , potentially because  often peaks late in the afternoon when VPD is highest.The   term was negligibly correlated with SWC in the riparian woodland (  = 0.02) but strongly correlated with SWC in the upland savanna (  = 0.76), likely  due to contrasting groundwater availability at the two sites (Mayes et al., 2020;Sabathier et al., 2021).The  term was the dominant driver of   at both sites (  = −0.73 and −0.66, respectively), which was expected given that  is encoded in the   calculations.SW↓ was the dominant driver of   (  = 0.86 and 0.91, respectively).The analysis of the environmental variables also explains the negative correlations between   −   and   at both sites (Fig. 7), which were contrary to expectations.The negative correlations likely emerge from the fact that  has negative correlations with both   and   (Table 1), yet   and   have opposing effects on   −   .Thus, the effect of  on   and   −   is likely large enough to confound the relationship between   and   −   .

Evaporative cooling
The modeled values of   and  , revealed the change in   that can be attributed to evaporative cooling of the canopy.The  , −   values were generally greater in the riparian woodland than the upland savanna (Fig. 8).At both sites, seasonal variability in  , −   tracked the seasonal trends of   .The smallest values of  , −   occurred in May at the riparian woodland and in May and June at the upland savanna.The largest values of  , −   (i.e., the most evaporative cooling) occurred in August at both sites.In the riparian woodland, the maximum daily climatological  , −  was 1.45 • C in May and 4.81 • C in August.In the upland savanna, the maximum daily climatological  , −   was 1.14 • C in May and 3.35 • C in August.The dip in  , −   values in the middle of the morning is likely a measurement or modeling artifact, potentially caused by shading of the flux tower sensors.

Impact of evaporative cooling on leaf respiration
Leaf light respiration ( ,ℎ ) was predicted using modeled values of   and  , .The difference between the two predictions (  ) indicates the change in  ,ℎ that is attributable to evaporative cooling of the canopy.In the riparian woodland,   exhibited consistent seasonal patterns each year, with the lowest values occurring during the pre-monsoon period in May and the largest values occurring in August (Fig. 9).In May, evaporative cooling decreased  ,ℎ by 5%-11%.In August, evaporative cooling decreased  ,ℎ by 21%-24%.In the upland savanna,   varied much more sporadically, likely due to the dependence of the ecosystem on water inputs from precipitation.The smallest values of   typically occurred in May and June of each year, and the largest values typically occurred in August.In May, evaporative cooling decreased  ,ℎ by 4%-11%.In August, evaporative cooling decreased  ,ℎ by 7%-28%.The largest value of   occurred in July 2008, when evaporative cooling decreased  ,ℎ by 31%.It is important to note the difference in sample size at the two study sites (3 years for the riparian woodland vs. 11 years for the upland savanna) due to the limited measurements of  ↓,  ↑, and  ↓ in the riparian woodland, which may account for some of the contrasting variability.

Discussion
We have presented a novel model to predict leaf temperature (  ) as a linearized function of the evaporative fraction (  ).The model predictions and empirical observations presented here demonstrate that evapotranspiration reduces   by consuming energy that would otherwise be partitioned into sensible heat flux.The model predicts that   −   varies as a linear function of   when all other variables are held constant.The model also predicts that   =   when   = 1.When   < 1,   theoretically varies as a function of   and   .The theoretical predictions from the energy balance model were tested using canopy-scale measurements of leaf temperature (  ) from two flux towers with contrasting water availability.At both sites,   converged to   when   approached 1.The mechanistic model presented in Eq. ( 7) exhibited strong model fit at both sites.Our findings are also consistent with a multi-site synthesis reported by Panwar et al. (2020), who demonstrated that the difference between surface temperature and air temperature was negatively correlated with   across a variety of ecosystems.

Environmental controls on 𝑇 𝑐
The flux tower observations suggest that water availability plays an important role in regulating   and its impact on   .The riparian woodland has consistent access to shallow groundwater (depth to groundwater ≈ 10 m), which provides a persistent source of water that is decoupled from the local precipitation regime on short time scales.The upland savanna does not have access to groundwater (depth to groundwater > 100 m) and is thus reliant on water inputs from monsoonal precipitation during the growing season.As a result,   was decoupled from near-surface SWC in the riparian woodland (  = 0.02) but strongly coupled to SWC in the upland savanna (  = 0.76).The enhanced water availability in the riparian woodland resulted in an earlier increase in   in the late spring and higher values of   throughout the growing season compared to the upland savanna (Fig. 6).
The differences in   at the two sites resulted in different magnitudes of evaporative cooling throughout the growing season.In the riparian woodland,   −  < 2 • C for much of the growing season, while in the upland savanna   consistently exceeded   by as much as 10 • C (Fig. 5).The seasonal patterns of   −   matched the seasonal patterns in   at both sites.That being said, the strength of the correlation between   −   and   in the riparian woodland decreased throughout the growing season, suggesting that   , and not water availability, was the primary driver of riparian   by the end of the growing season (Fig. 7).In the upland savanna,   −   was most strongly correlated with   in the middle of the growing season during peak monsoonal precipitation and less strongly correlated with   during the drier periods at the beginning and end of the growing season.
At low values of   ,   −   is largely regulated by non-evaporative cooling processes (Muller et al., 2021(Muller et al., , 2023)).The efficiency of nonevaporative cooling is determined by the resistance to sensible heat flux (  ).At a leaf scale,   is a function of leaf size, leaf structure, and the wind speed across the leaf surface (Balding and Cunningham, 1976;Jones, 2014;Leigh et al., 2017).At a canopy scale,   is also a function of vegetation cover and vegetation structure, which drive turbulent mixing (Yang and Friedl, 2003;Rigden et al., 2018).The upland savanna experienced a larger range of   −  values at low values of   because there was a larger range of   values under those conditions.In the riparian woodland,   and   covaried more strongly, resulting in a smaller range of   −   values at low values of   .Interestingly, the upland savanna experienced higher wind speeds and likely had smaller leaves (Stromberg et al., 1993), which are typically associated with more efficient heat transfer, but the riparian woodland had lower modeled values of   .This suggests that turbulent mixing at a canopy scale played an important role in regulating   , which is consistent with previous analyses of   across different vegetation types (Rigden et al., 2018;Young et al., 2021).

Limited homeothermy
The mechanistic relationship between   and   has received considerable attention in literature (Cavaleri, 2020, and references therein), with various studies arguing that plants exhibit either limited homeothermy (  <   at high values of   ), poikilothermy (  ≅   ), or megathermy (  >   at high values of   ).We found that the riparian woodland generally exhibited poikilothermy.The slope of the relationship between   and   was close to 1 ( = 0.92) and there were few observations where   <   , even at high values of   .The upland savanna exhibited megathermy; the slope of the relationship between   and   was greater than 1 ( = 1.17) and   consistently exceeded   .Neither site in this study exhibited a clear signal of limited homeothermy.Moreover, the mechanistic model of   always predicts that   ≥   when   ≥ 0 W m -2 and   ≤ 1.Even if stomatal conductance is not limiting, there is by definition not enough   in the system to increase  to levels that result in   <   under normal conditions.It follows from Eq. ( 6) that   <   can only occur if   > 1.Previous studies have demonstrated that   > 1 only occurs briefly around sunrise and sunset when  is negative and  is positive, a time of day when the magnitudes of energy fluxes are small.The value of   is somewhat constant during daylight hours and typically substantially less than 1 (Crago, 1996;Gentine et al., 2007Gentine et al., , 2011)).Conditions where   > 1 can also occur as a result of the ''oasis effect'' whereby the advection of dry air over well-watered vegetation creates a land-atmosphere feedback that causes  to exceed   (Baldocchi et al., 2016).The oasis effect is most commonly associated with rice paddies and wetlands in semi-arid climates, but it is not clear how often the effect actually occurs (Baldocchi et al., 2016).
Despite the lack of theoretical or empirical support for limited homeothermy in the data examined here, observations where   <   are commonly reported in literature.Some researchers have suggested that observations of limited homeothermy are due to systematic errors from certain types of in situ sensors (Still et al., 2019b).However, observations where   <   have also been reported in studies that measure   using infrared radiometers (e.g., Idso et al., 1981;Jackson et al., 1981;Kar and Kumar, 2007;Ballester et al., 2013;Blonder et al., 2020).Thus, there is an apparent paradox whereby observations of   <   seem highly unlikely given fundamental energy balance constraints (Eq.( 6)), but are nonetheless common.Blonder and Michaletz (2018) demonstrated from energy balance theory that limited homeothermy can only occur when stomatal conductance is high and   is low.Other research has examined non-steady state   dynamics, which are not explored here (e.g., Leigh et al., 2017).The relationship between   and   established by this study represents another novel constraint on leaf thermoregulation via limited homeothermy.Further theoretical and empirical research is needed to constrain the conditions that result in observations where   <   , especially given the substantial disagreement over how frequently leaf thermoregulation actually occurs in nature (e.g., Blonder et al., 2020;Still et al., 2022).

Plant carbon balance
Constraining the mechanistic relationship between   and   is of critical importance for modeling ecosystem responses to anthropogenic climate change.Leaf energy balance and   serve as fundamental constraints on the selection and adaptation of plant traits (Michaletz et al., 2015(Michaletz et al., , 2016)), which are generally assumed to maximize net carbon uptake while controlling for the risk of plant hydraulic failure (Wolf et al., 2016;Sperry et al., 2017;Mencuccini et al., 2019).Previous trait-based research has focused on the role of stomatal conductance in maximizing photosynthetic assimilation via biochemical fixation of carbon (Cowan and Farquhar, 1977;Medlyn et al., 2011).Here we demonstrate that stomatal conductance also alters net carbon uptake via the impact of evaporative cooling on   .In both the riparian woodland and upland savanna, evaporative cooling of the leaf surface often reduced   by ca.15% in the middle of the growing season.Reduced   from evaporative cooling would also be expected to keep   closer to the photosynthetic optimum in hot environments, maximizing photosynthetic assimilation (Roden and Pearcy, 1993;Medlyn et al., 2002).For example, Uni et al. (2022) demonstrated that a reduction in   from 40 • C to 35 • C would increase photosynthetic assimilation by 42%.Their study analyzed Acacia tortilis, a species that is structurally and functionally similar to P. velutina.The regulation of   by stomatal conductance represents an important linkage between plant water and carbon cycles that has received little attention in literature (but see Michaletz et al., 2015Michaletz et al., , 2016) ) and may alter predictions of optimal plant traits and behavior.All other factors held constant, the data examined here suggest that high levels of  will enhance net carbon uptake by reducing   , which may marginally favor high risk-high reward hydraulic strategies in dryland vegetation (e.g., Hultine et al., 2020;Williams et al., 2022).

Thermal remote sensing
Eq. ( 7) also provides a physical basis to interpret thermal remote sensing measurements (Mallick et al., 2022).Tower-mounted infrared radiometers are a reliable proxy for airborne and satellite thermal sensors, which can measure surface temperature over broad spatial scales.Thermal remote sensing is widely used to monitor agricultural productivity (Jones et al., 2009;Maes and Steppe, 2012) and manage water resources (Anderson et al., 2012).Our study joins other recent efforts to unify plant traits and thermal measurements, which will likely yield novel insights into ecosystem processes at leaf to global scales (Still et al., 2019a(Still et al., , 2021;;Farella et al., 2022).

Conclusion
The mechanistic relationships between water, energy, and carbon fluxes at the leaf surface are of considerable importance for predicting the responses of terrestrial ecosystems to anthropogenic climate change.The model presented here constrains the mechanistic relationship between   and   and provides a framework to quantify evaporative cooling of the leaf surface.Importantly, the model reveals that   −  varies as a linear function of   and that   −  = 0 • C when   = 1.The model predictions were validated using measurements of canopy-scale leaf temperature (  ) from two flux towers.Seasonal variability in measured   was primarily driven by   , although   also played an important role in regulating   in well-watered conditions.Neither the model predictions nor the empirical observations provided evidence for regimes where   is substantially less than   .Future work is needed to understand the conditions that result in empirical observations of   <   in croplands.Our analysis also reveals that evaporative cooling of the leaf surface has important functional implications for plant carbon cycling.Evaporative cooling substantially reduced   at both study sites.The impact of evaporative cooling on   may affect predictions of optimal plant traits and behavior under future climate scenarios.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 3 .
Fig. 3. Boxplots of daily  ∕ values in the riparian woodland (a) and upland savanna (b) for each month of the growing season.Outliers are not shown.Source: Data are reanalyzed from Scott et al. (2021) and Nelson et al. (2020a,b).

Fig. 4 .
Fig. 4. Energy flux measurements from the riparian woodland (a,c) and upland savanna (b,d) differ in that the riparian woodland has a higher evaporative fraction (  ) compared to the upland savanna.However, in both systems leaf temperature (  ) converges to air temperature (  ) when   approaches 1.The plots in the top row (a,b) compare   and   measurements.The blue lines are the 1:1 line.The plots in the bottom row (b,d) show the relationship between   −   and   .White areas in the plots indicate that there were 0 observations in that portion of the feature space.The color scale saturates when there are more than 75 half-hourly observations in a given portion of the feature space.Some outliers are outside of the plotted range.

Fig.
Fig. Diurnal and seasonal climatology of   −   (a,b) and   (c,d) in the riparian woodland (a,c) and upland savanna (b,d).The colored lines represent the mean values for each time of day, grouped by month of the growing season.A comparable figure for   ,   , and   is located in the Supplementary Materials.

Fig. 6 .
Fig. 6.Seasonal climatology of evaporative fraction (  ), resistance to sensible heat flux (  ), and available energy (  ) in the riparian woodland (a-c) and upland savanna (d-f).The colored bars represent the median values for each month of the growing season.The black vertical bars represent the interquartile range.The   values are limited to observations from 12:00 local time.

Fig. 7 .
Fig. 7. Spearman rank correlations between   −   and evaporative fraction (  ), resistance to sensible heat flux (  ), and available energy (  ) at the riparian woodland (a-e) and upland savanna (f-j).The colored bars represent the correlations for each month of the growing season.

Fig. 8 .
Fig. 8.Diurnal and seasonal climatology of modeled  , −   , which indicates the change in   due to evaporative cooling of the leaf surface.The colored lines represent the mean values for each time of day, grouped by month of the growing season, for the riparian woodland (a) and upland savanna (b).

Fig. 9 .
Fig. 9. Monthly mean decrease in daytime leaf respiration (  ) that is attributable to evaporative cooling for the riparian woodland (a) and upland savanna (b).The black vertical bars indicate the range of monthly mean   values for individual years.There are 3 years of data for the riparian woodland and 11 years of data for the upland savanna.

Table 1
Spearman rank correlations between observed   ,   , and   and environmental variables measured by the flux towers, including shortwave insolation (SW↓), vapor pressure deficit (VPD), wind speed (), and soil water content (SWC) for the riparian woodland and upland savanna.