Studying the large-scale effect of leaf thermoregulation using an Earth system model

Plants have the ability to regulate heat and water losses. This process also known as “leaf thermoregulation” helps to maintain the leaf temperature within an optimal range. In a number of laboratory and field experiments, the leaf temperature has been found to deviate substantially from the ambient temperature. In the present study, we address the question of whether the negative correlation between the leaf temperature excess and the ambient air temperature, which is characteristic of leaf thermoregulation, constitutes a robust feature at larger scales, across a broad range of atmospheric conditions and canopy 5 characteristics. To this end, we developed a new dual-source canopy layer energy balance scheme (CEBa) and implemented it into JSBACH, the land component of the Max Planck Institute for Meteorology’s Earth system model (MPI-ESM). The approach calculates the temperature and humidity in the ambient canopy air space, the temperature of the ground surface, and the temperature of the leaf as well as the energy and moisture fluxes between the different compartments. Here leaf thermoregulation is investigated using different modeling approaches, namely a zero-dimensional instantaneous solution of the 10 energy balance as well as offline FLUXNET site experiments and coupled global simulations. With the help of the simulations at the site-level, we can show that the model is capable of reproducing the effect of leaf thermoregulation even though the simulated signal at the canopy scale is less pronounced than indicated by measurements at the leaf scale. However, on a global scale and over longer-timescales, this negative correlation is only simulated in idealized setups that neglect limitations on the plant available water, and even then, the signal is less pronounced than indicated by the short-term observations of individual 15 leaves. When accounting for moisture limitations, we predominantly find positive correlations between leaf temperature excess and the ambient air temperature.


Introduction
The canopy is commonly defined as "the community of aboveground plant organs" (Campbell et al., 1989) with over two-20 thirds of the terrestrial surface being covered by vegetation (Huete et al., 2004). The canopy influences the turbulent exchange of energy, water, and momentum between the land surface and the overlying air mass via a number of biogeophysical and biogeochemical processes. Furthermore, its complex structure and the resultant shading strongly affect the terrestrial radiation budget, largely regulating the distribution of incoming solar energy into ground and turbulent heat fluxes. This gives the canopy a key role in controlling surface and below-ground temperatures, which in turn affect the diurnal variation of the boundary layer 25 development and important atmospheric processes such as cloud formation and convection.
Moreover, the canopy determines key characteristics of the hydrological cycle, via a strong influence on evaporation and transpiration (e.g., Oki and Kanae, 2006). This affects dynamics and thermodynamics of the climate system (Chahine, 1992), especially because soil moisture is part of several feedbacks on the local, regional, and even on the global scale (Seneviratne et al., 2010;Lawrence et al., 2007). The canopy also constitutes an important element of the global carbon cycle, as substantial 30 amounts of atmospheric carbon are taken up by the leaves within the canopy layer. Woody plants alone assimilate about 120 Gt carbon per year through the process of photosynthesis (e.g., Schimel et al., 2001), with net primary productivity (NPP) being largely constrained by atmospheric CO 2 concentrations and the available soil moisture, but also by the prevailing leaf temperatures (e.g., Farquhar et al., 1980;Farquhar and Sharkey, 1982;Ball et al., 1987).
The first detailed studies on the energy balance of leaves and the heat transfer between the leaf and its environment date 35 back to the 1960s. Here Gates (1962) was amongst the first to show that leaf temperatures can deviate substantially from the temperature of the surrounding air: Sunlit leaves were up to 20 • C warmer than the ambient air, while shaded leaves were found to be on average 1.5 • C colder than the air temperature. These temperature differences occur on time scales ranging from seconds to minutes and can vary significantly throughout the day due to varying environmental conditions such as solar and thermal radiation, air temperature, wind speed, and water vapor content (Gates, 1965). Linacre (1964) introduced the concept 40 of "leaf thermoregulation" (LT), demonstrating a distinct negative correlation between the ambient air temperature and the leaf temperature excess -which is the difference between leaf and ambient air temperature -for water-unstressed, photosynthetic plants ( Fig. 1). This study also showed that the leaf temperature excess changes its sign, allowing the leaf to be warmer than the surrounding air in cold environments, and to stay colder than the air in warm environments. The temperature at which this occurs is called "equivalence temperature", which Linacre (1967) measured to be around 34 • C. Linacre hypothesized that the

CEBa
The standard version of JSBACH (hereafter called Classic) and the recently implemented, extended SkIn + scheme (Heidkamp et al., 2018) both employ a closure of the surface energy balance, where the vegetation is represented by a single layer, based on the big-leaf approach. In these approaches, only one temperature -the surface temperature, corresponding roughly to the temperature at the displacement height -is used to characterize the canopy air space, the ground under or next to the canopy, 100 and the leaves within the canopy. To be able to represent leaf thermoregulation (LT) with JSBACH, we implemented the CEBa scheme which allows the calculation of four prognostic variables: the leaf temperature T leaf , the ground temperature T grd , the temperature in the canopy air space (CAS) T cas and the specific humidity in the CAS q cas . The canopy temperature T c is often described by the temperature of the entire above-ground vegetation, including stems and branches. Since CEBa does not distinguish between different parts of the biomass in the canopy, in particular the thermodynamic properties of branches, 105 trunks, and leaves, which are all described by the same temperature, we use T c as a proxy for the leaf temperature T leaf .
Many land surface schemes with a dual-source canopy layer (e.g., Samuelsson et al., 2006) solve the system of energy balance equations diagnostically for T cas and q cas without accounting for the heat capacity of the canopy air space, neglecting thermal energy storage as well as changes in enthalpy due to changes in the composition of the CAS. In contrast, we follow Vidale and Stöckli (2005) and consider heat and water storage in the CAS, which can not be neglected in the case of tall 110 vegetation (Heidkamp et al., 2018). In CEBa, the biomass heat storage depends on the vegetation temperature, while T cas is used to calculate the heat storage of the moist air in the canopy. Here the latent heat storage (i.e., heat contained in the humidity in the canopy layer) is no longer calculated by means of the effective surface specific humidity (Heidkamp et al., 2018). It is now calculated more realistically using the specific humidity q cas , representing the moisture content in the CAS.
In Classic and SkIn + , the distribution of the vegetation follows the so-called "land mosaic" configuration, in which the 115 plants are assumed to be clumped in one subgrid area of a grid box (vegetation fraction) and the bare soil occupies the rest (non-vegetation fraction). Both fractions receive the same amount of radiation and the total evaporation is calculated using area-weighted means. In contrast, in CEBa a sparse canopy is used in which the vegetation is distributed uniformly, similar to a savanna-type landscape (Lee, 2018). Here the total evaporation is the sum of transpiration and bare soil evaporation, while the weighting of their distribution is achieved indirectly by partitioning the net radiation R net using the so-called "sky-view 120 factor" χ (Verseghy et al., 1993) Here f veg is the vegetated fraction of the grid box and LAI eff is the effective leaf area index taking into account the clumping of vegetation by canopy gaps (Loew et al., 2014). This leads to a radiation flux of the canopy R net,c and of the underlying ground R net,g (Fig. 2). Note that in the following equations the index c pertains to the canopy (which serves in this study analogously 125 for the leaf) and the index g to the ground beneath: Here α is the albedo, ε the surface emissivity, σ the Stefan-Boltzmann constant, L in the incoming longwave radiation and S in the incoming shortwave radiation. The factor 2 in R net,c follows the assumption that the canopy layer emits radiation in two different directions: to the sky and to the soil surface. The radiation budget leads to two different cases of limit value 130 calculation: For a negligible leaf area index (LAI = 0) the sky view factor is maximal (χ = 1), which means that there would be no canopy at all leading to the known solution for bare soil: In contrast, a maximum leaf area index (LAI → ∞) results in a minimum sky view factor of χ = 0 and an opaque canopy layer that absorbs all incoming shortwave radiation: Here the soil surface still absorbs and emits longwave radiation in interaction with the overlying dense canopy layer. Figure 2 illustrates the energy and water exchange between the canopy, the ground, the canopy air space, and the lowest atmospheric level (LAL) of the new CEBa system, including the temperatures, humidities, and resistances, respectively. The sensible and latent heat fluxes are parameterized using the common resistance formulation.
(1-χ)εσT leaf 4 (1-χ)εσT grd 4 χεσT grd 4 RH up Figure 2. Principal sketch of the dual-source canopy energy balance scheme CEBa: Depicted are the radiative heat exchanges (yellow for shortwave radiation, orange for longwave radiation), temperatures (red), specific humidities (blue), energy fluxes (red arrows), and water fluxes (blue arrows), including the resistances (black), where χ is the sky-view factor, Lin the incoming longwave radiation, Snet the net shortwave radiation, ε the surface emissivity, σ the Stefan-Boltzmann constant, Tair the air temperature, T leaf the leaf temperature, T grd the ground surface temperature, Tcas and qcas the temperature and specific humidity in the canopy air space, qair the specific humidity of the air, qsat the saturated specific humidity at the respective temperature, ratm the atmospheric resistance, r b the bulk resistance, r d the aerodynamic resistance between ground and CAS, rc the canopy resistance, and RHup a nonlinear function depending on the relative humidity of the uppermost soil layer.
There are three different pairs of sensible and latent heat fluxes that are connected to each other through the CAS. The first is the exchange of energy and water between the CAS and the atmosphere: where ρ is the density of humid air, c p the specific heat capacity of air, L v the specific enthalpy of vaporization, T air and q air the temperature and specific humidity of the air at the LAL, T cas and q cas the temperature and specific humidity of the air in 145 the CAS, and r atm the atmospheric aerodynamic resistance, which is the inverse product of wind speed and drag coefficient.
The drag coefficient represents a measure of the turbulence strength depending on the roughness of the underlying surface and the atmospheric stratification. The second pair is the sensible and latent heat flux between the CAS and the canopy biomass (for simplicity's sake snow and interception evaporation is not considered here): where T c is the temperature of the canopy surface (here a proxy for the leaf temperature T leaf ), q sat the saturated specific humidity at the respective temperature, r b the aerodynamic resistance between the canopy and the CAS, and r c the canopy resistance. The latter describes the additional resistance that water vapor molecules encounter when moving through the stomatal openings of leaves. It is the reciprocal of the stomatal conductance, which depends on the assimilation rate and hence on the carbon flux through the stomata (Farquhar and Sharkey, 1982). It is scaled by the LAI (leaf area index) to achieve the transition 155 from leaf scale to canopy scale. In addition, it is modified by a water stress factor that depends on the saturation of the root zone. Since photosynthesis depends on the leaf temperature, in CEBa, r c depends on the leaf temperature instead of the air temperature of the LAL. The resistance between the leaf and the canopy air space r b (also called bulk resistance) is a function of the LAI, the leaf width, and the wind speed at the top of the canopy u cas and describes the additional resistance induced by the turbulence in the leaf boundary layer (Choudhury and Monteith, 1988). It is modified with a free convection correction 160 according to Sellers et al. (1986) depending on T c , T cas and the LAI. The last pair is the sensible and latent heat flux between CAS and the ground surface under the canopy: where T g is the temperature of the ground surface, r d the aerodynamic resistance between the ground surface and the CAS, and RH up a nonlinear function depending on the relative humidity of the uppermost soil layer. The latter serves as a resistance 165 for bare soil evaporation. The parameterization of r d follows a more complex approach depending on the vegetation height z veg , the zero plane displacement height d (again a function of z veg ), and the wind speed at the top of the vegetation u cas . It is based on the general equation of an aerodynamic resistance for momentum transfer between the soil surface and the sink for momentum in the vegetation under neutral conditions (Choudhury and Monteith, 1988). We also apply a stability correction for unstable cases following Sellers et al. (1986). For the estimation of u cas , we use an approach according to Xue et al. (1991).

170
The basic idea is the concept of a roughness layer defined by a so-called transition height z trans (see Appendix B).
Considering these expressions, we obtain a system of four coupled equations: Here C c and C T are the area-specific heat capacities of the canopy and the CAS, respectively, L q is the area-specific enthalpy of vaporization for water, G is the ground heat flux, F CO2 is the net CO 2 flux and β = 10.884 · 10 6 J/kg is a conversion factor 175 corresponding to the energy required to incorporate 1 mol CO 2 by carbohydrate bonds through the process of photosynthesis.
We assume an infinitesimally thin surface at the ground, which allows to include the ground heat storage by the ground heat flux G. The nonlinear terms within the system of equations in CEBa are linearized using a first-order Taylor approximation and the time derivatives are discretized. The resulting linear system of four prognostic variables is solved by matrix inversion and iteration.

Experiments and data
For the first experiment, we do not run the full ESM or even the complete land-surface component. Rather, we solve the system of energy balance equations of the CEBa model for different air temperatures ranging from −20 • C to 60 • C. Here we assume an unlimited soil water availability for transpiration (no water stress), as well as for bare soil evaporation (RH up = 1). The surface parameters were chosen to represent a forest with tall vegetation -i.e., z veg = 30 m, z 0 = 4.8 m, LAI = 4.5 m 2 /m 2 , 185 where z veg is the vegetation height, z 0 the roughness length, and LAI the leaf area index. Furthermore, we prescribe the specific humidity of the air, as well as the wind speed and the incoming solar radiation. By varying the latter, we can analyze how the equivalence temperature (T eq ) and the slope of the regression line (SRL) depend on the relative humidity, incoming solar radiation, canopy resistance, and aerodynamic leaf resistance. Although the leaf temperature excess as a function of air temperature represents a nonlinear curve, we approximate it by a linear function, in order to describe the effect of LT simply 190 by the SRL and T eq .
In this experiment, we calculate the stationary solution as a function of the air temperature to determine T eq and the SRL describing the negative correlation due to LT. This steady-state experiment can be interpreted as a zero-dimensional experiment since the solution neither depends on space nor on time and there are neither heat nor water mass storages. In this case, the heat fluxes between the CAS and the LAL are the sums of the ground and canopy heat fluxes. Thus, the third and fourth equation 195 of (8) simplify to 8 https://doi.org/10.5194/esd-2020-75 Preprint. Discussion started: 19 November 2020 c Author(s) 2020. CC BY 4.0 License.  To obtain the total net radiation the partial radiative fluxes have to be weighted by the sky view factor χ: Note that it is debatable whether air temperature measurements resemble T cas , the temperature next to (or even under) the 200 vegetation within the CAS, or T air , the air temperature above the canopy. In the respective analysis, the leaf temperature excess is defined as the temperature difference between T leaf and T air .
To address the second scientific question, i.e. whether a negative correlation of leaf temperature excess and air temperature emerges under realistic conditions at the canopy scale, we performed two offline single-site experiments with JSBACH including CEBa, using 30-minute measurement data from a FLUXNET tower in Brazil's tropical forest (Goulden, 2016) and 205 the temperate Tharandt forest in Germany (Bernhofer et al., 2016). Both datasets are part of the "FLUXNET2015 Dataset" (Pastorello et al., 2020). In an offline experiment, the LSM is decoupled from its host model and forced by observation data (precipitation, air pressure, air temperature, specific humidity, and wind, as well as shortwave and longwave incoming radiation).
Both are forest sites with a dense canopy and a large vegetation height. The tropical site (ID: To investigate the effect of LT on the global scale, we performed a set of coupled simulations with the MPI-ESM at T63 resolution (i.e. 1.9 • ), covering the years 1979 to 2008. The simulations follow the AMIP protocol (Gates, 1992), with the sea surface temperatures and the sea-ice extent being prescribed. Technically, it is possible to run the CEBa scheme for shallow or even no vegetation, such as deserts or glaciers (the latter with an average vegetation height of z veg = 0). However, the 225 simulations with CEBa for these extremes produce implausible results. Therefore, the CEBa scheme is used only in grid boxes with an average vegetation height greater than 5 m. The resulting mask is depicted in Fig. 5. For all other land points, the SkIn + scheme is used and the grid boxes are excluded from the following analysis. Both the offline and coupled simulations are run with a time step of 7.5 min.

230
In this section the process of leaf thermoregulation (LT) is analyzed as follows: First, we asses the influence of atmospheric conditions and leaf properties on the leaf's ability to regulate its temperature. Secondly, we perform offline experiments with JSBACH including the CEBa scheme for a tropical and a temperate forest site to identify the effect of LT at the canopy scale under realistic conditions. Finally, we perform a global coupled model experiment with CEBa to establish whether a negative correlation between leaf temperature excess and ambient air temperature can be found over a sizeable spatial area and on a 235 decadal time scale in model simulations.

Steady state experiment
In the temperature range between 0 • C and 55 • C, the system of energy balance equations of the new CEBa scheme (Eq. 8) reproduces the negative correlation between leaf temperature excess ∆T and ambient air temperature that is characteristic of LT (Fig. 6, see also Appendix C). For low air temperatures, where transpiration is negligible, the leaf is warmer than the 240 surrounding air because the incoming shortwave and longwave radiation is absorbed and transformed into internal energy, resulting in an increase of the leaf surface temperature. In equilibrium, this energy is released in the form of fluxes of outgoing longwave radiation and sensible heat, which transfer energy from the leaf surface to the air and compensate for the radiative forcing. For higher temperatures, ∆T decreases almost linearly with increasing air temperatures, which is a result of the nonlinear temperature-dependency of the latent heat flux. The sensible and latent heat fluxes are mainly driven by the temperature 245 and moisture gradient between leaf and air and for a linear temperature increase, the saturated specific humidity increases exponentially following Clausius-Clapeyron theory. At a (comparatively) constant net radiation and relative humidity in the canopy and assuming water-unstressed plants, this results in an energetic equilibrium, in which the latent heat flux increases with rising temperatures, whereas the sensible heat flux decreases. Above a certain temperature, the latent heat flux exceeds The slope of the regression line (SRL) and T eq -which can be used to characterize the effect of LT -depend on atmospheric conditions, such as incoming radiation, but also on the characteristics of the canopy, such as the canopy resistance (Fig. 7).

255
T eq increases for higher radiative fluxes and for higher humidity, with the T eq dependency on the incoming radiation following a logarithmic function and an exponential function for relative humidity. This implies that variations in solar radiation are especially important when the solar incoming radiation is low (lower than 400 W/m 2 ), while variations in the air's moisture content are important when the relative humidity is high (higher than 80 %). The SRL increases with saturation but decreases with rising radiative fluxes, indicating LT to be more effective at higher radiative fluxes and in dryer atmospheric conditions.

260
Here the SRL displays an almost linear dependency on incoming solar radiation and relative humidity, indicating a similar sensitivity at high and at low values. In contrast, changes in canopy resistance r c and aerodynamic leaf resistance r b exhibit a notable influence on the LT curve almost exclusively at low values of r b and r c . Here T eq increases logarithmically with increasing resistances for both r c and r b . The SRL is substantially affected by changes in r c and r b when assuming a high incoming solar radiation, while for lower radiative fluxes the SRL is largely independent of the resistances.   Michaletz et al. (2016). The brightness of the colors indicates how many percents of the data fall within an "area" of 1 K × 1 K = 1 K 2 .

FLUXNET Sites
The above results show that the strength of the LT effect -given by the SRL -is in large parts determined by the atmospheric conditions, with a strong LT being possible only in case of a dry atmosphere and a strong solar flux. This raises the question whether the negative correlation between leaf temperature excess and air temperature also constitutes a robust feature under realistic atmospheric conditions. To adress this question, we use JSBACH with the model being forced using observations 270 from two FLUXNET sites located in the temperate needleleaf Tharandt forest in Germany and in a tropical evergreen forest in Brazil.
At the Tharandt site, the daytime air temperatures during the vegetation period range between 0 • C and 35 • C and follow a Gaussian distribution (Fig. 8a). The leaf temperature excess shows a sharp drop in frequency for negative temperature differences and the two-dimensional probability function indicates a larger spread towards positive (3.7 K for frequencies larger than 275 0.25 %/K 2 ) than towards negative leaf temperature excesses (−1.6 K for frequencies larger than 0.25 %/K 2 ). Nevertheless, the regression line between leaf temperature excess and air temperature exhibits a negative correlation, showing that the effect of LT is present also at the canopy scale under observed atmospheric conditions. Here T eq is equal to 31.7 • C and similar to the 30.1 • C reported by Michaletz et al. (2016). However, the SRL (−0.1) is significantly smaller than the SRL of −0.27 found by Michaletz et al. (2016) and the simulated LT is constantly weaker than indicated by observations. One potential reason for 280 the mismatch between the simulated and the observed SRL is the water-stress level of the vegetation. While the simulations account for limitations due to water-availability, Michaletz et al. (2016) investigated water-unstressed plants. When neglecting water stress in the simulations, the distribution of ∆T shifts towards negative temperature differences -the largest significant (> 0.25 %/K 2 ) negative differences now amount to −3 K -and the simulated SRL (−0.24) agrees well with the observations (Fig. 8b). However, the stronger evaporative cooling reduces the simulated T eq to 18.6 • C, which is around 12 K below the 285 value reported by Michaletz et al. (2016).
In general, the simulated leaf temperatures are substantially lower than those reported by Michaletz et al. (2016) resulting in an almost constant offset between the simulated and observed regression lines in the case of unstressed vegetation (Fig. 8b). This is mainly due to focus and scale differences between the observations and the simulations. Most temperature measurements target individual leaves at the surface of the canopy, while the CEBa scheme covers the entire extent of the canopy layer.

290
Thus, while the observations predominantly pertain to the temperatures of the warmer leaves, which receive direct sunlight, the CEBa scheme estimates the average over the entire canopy biomass, including the shaded parts of the canopy, which are substantially cooler.
At the tropical site, the simulated effect of LT is not only weaker than the observed values, but the leaf temperature excess even shows a positive correlation (SRL = 0.15) with the air temperatures (Fig. 9a). Similar to the Tharandt site, one reason for 295 the high SRL is the increasing water stress at higher temperatures, and neglecting the limitations due to soil water availability results in a negative correlation between ∆T and T air (Fig. 9b). However, even under unstressed conditions, the SRL (−0.18) is substantially higher than the observed values (−0.27) or the SRL simulated for the Tharandt site (−0.24), indicating a weaker LT effect in tropical forests. Here the reason for the weak LT lies in the prevailing high relative humidity in the atmosphere, which limits the potential for an evaporative cooling of the leaf. While in the unstressed case the average relative humidity 300 at the Tharandt site was 62 %, the tropical site featured a relative humidity of 71 %. Thus, while we find that the strength of LT at the canopy scale under realistic atmospheric conditions is comparable to observations at the leaf scale, this is only the case for dry atmospheric conditions and under water-unstressed conditions. Once limitations due to soil water availability are accounted for, the negative correlation between ∆T and T air becomes weaker and even becomes positive for high levels of relative humidity.

Global AMIP experiments
The positive correlation between ∆T and T air is not unique to the investigated site in the tropics but we also find it to be the dominant long-term signal at the global scale. In an AMIP experiment, the 30-year mean daytime leaf temperature during the vegetation period (Apr-Sep) is consistently larger than the corresponding air temperature, with the smallest ∆T amounting to 0.5 K. Here we assume that a space-for-time-substitution is justified. This means that the variations of the data points (30-year 310 averages) in Fig. 10, which result due to variations in space at different grid points, can be compared with those data points measured at a single site that are a function of time at a fixed location. Large ∆T are predominantly simulated in grid cells that also feature high air temperatures, which leads to a positive correlation between leaf temperature excess and air temperature, with a SRL of about 0.04 (Fig. 10, red circles). to water-availability, where a larger water stress is predominantly simulated for grid cells with higher air temperatures. An idealized experiment without water-stress (where unlimited soil water availability is assumed) produces a negative correlation between ∆T and T air , with a SRL of around -0.01 (Fig. 10, blue circles). Here ∆T constitutes predominantly positive values, which is in good agreement with the observations. However, these positive ∆T are significantly smaller than the short-term measurements. Negative temperature difference are hardly present, but a slight negative trend is visible at air temperatures of 320 about 30 • C, which would agree with the T eq reported by Michaletz et al. (2016). However, the model does not even simulate 30-year means of daytime air temperature higher than 30 • C. This leads to a significantly smaller SRL than the values derived in the site-level experiments of about -0.01. One reason for this feature could be caused by a negative land-atmosphere-feedback, which is present in the AMIP experiment, but not in the site level simulations, where the air temperature and specific humidity at the lowest atmospheric model level are prescribed based on observations. In the AMIP experiment, the evaporative cooling 325 does not only reduce the leaf temperature but also affects the atmospheric temperatures, diminishing the overall effect on ∆T .
The theory of LT postulates that leaves adjust their temperature to maximize the carbon gain of the plant, but the respective increase in productivity has not yet been assessed at the global scale. In the following we compare the gross primary productivity (GPP) of the above AMIP simulations with that of an additional set of simulations in which the GPP and the canopy resistance are simulated based on air temperatures. Here we analyze the correlation between the GPP (in µmol CO 2 /(m 2 s)) calculated 330 based on air temperatures GPP(T air ) and the difference ∆GPP = GPP(T leaf ) − GPP(T air ) for the water-stressed and waterunstressed case, with positive values of ∆GPP indicating an increase in GPP due to the consideration of LT. Both the stressed and unstressed experiment show a general increase in GPP with regression curves above the x-axis (Fig. 11). However, while the unstressed experiment shows almost exclusively positive values, the stressed experiment simulates a negative ∆T in several cases. In general the simulated GPP differences are rather small and the effect of LT appears to be strongest for GPP rates of 335 around 4 to 8 µmol CO 2 /(m 2 s) with a maximum of the regression curve of 0.28 µmol CO 2 /(m 2 s) at 5.1 µmol CO 2 /(m 2 s).
Overall, using the leaf temperature instead of the air temperature to calculate the canopy resistance and GPP leads to an increase in the 30-year global GPP mean of around 3 % for the water-stressed experiment and 5 % for the unstressed case. Thus, our results support the theory that the temperature regulation of the leaf constitutes a strategy to maximize the carbon gain of the plant.

4 Discussion and Conclusion
This study constitutes the first investigation of the effect of leaf thermoregulation (LT) at larger spatial and temporal scales.
We demonstrated that the characteristic correlation between leaf temperature excess and air temperature is not only a feature of individual leaves but can also be found at the canopy scale. Here our model simulated a negative correlation on the site- Our results indicate that the leaf temperature excess is substantially affected by the prevailing soil moisture conditions. The

355
representation of water stress in JSBACH has not yet been evaluated explicitly and could still contain uncertainties, e.g., due to biases in the soil moisture content. Moreover, there are potential water limitations that could control the canopy resistance, possibly, due to missing processes. The temperature dependence for the calculation of the stomatal conductance as well as for the photosynthesis rate is assumed to follow an Arrhenius-type temperature dependence. This means in the context of JSBACH: the higher the temperature, the higher the stomatal conductance or the larger the GPP rate (up to a certain temperature where a 360 high-temperature inhibition sets it to zero). What is missing so far is a parameterization of a process that artificially regulates the stomata to the plant's temperature optimum.
We have shown that the effect of LT (i.e., the amplitude of the leaf temperature excess) depends to a large extent on the relative humidity of the ambient air and the increasing solar insolation. Regarding the latter, it is crucial to estimate the distribution of radiation within the canopy layer correctly, i.e., how much radiation is absorbed or reflected from the leaf and how much 365 is transmitted to the ground surface. In this context, studies have shown that leaves have the ability to reduce or enhance the absorption of solar radiation, e.g., by increasing the leaf or needle density, by rotating the leaf angle, or by reflective leaf hairs (Helliker and Richter, 2008).
We found that, especially during midday, the leaf temperature is almost always higher than the air temperature since it is completely illuminated. In this case, if the leaf temperature were lower than the air temperature at high solar radiation, 370 this would represent a state with stable stratification, which is not realistic in convective conditions during the day. To avoid that, only a large heat flux from the ground surface could compensate that imbalance. Thus, it is questionable whether LT is even possible for extremely dense or closed canopy layers (e.g., in the tropics). In that regard, the implementation of a flux aggregation method to solve the energy balance equation (Best et al., 2004;de Vrese and Hagemann, 2016) could provide a more realistic representation of the canopy layer.

375
Our results demonstrate that future studies need to consider the vertical distribution of leaf temperatures regarding LT theory to distinguish correctly between leaf and canopy scale. In this context, a multilayer model is probably more appropriate as it can use parameters measured at the leaf scale and at different heights to calculate the vertical temperature gradient in the canopy layer (Dai et al., 2004). Luo et al. (2018) suggest to avoid the big-leaf approach and advocates the use of a two-leaf from the leaf to the canopy scale.
For exploring the role of the leaf temperature in the canopy, we made a first approach to investigate the effect of LT at the canopy scale using an ESM. We show that we are indeed able to detect a negative correlation between leaf temperature excess and air temperature, but its magnitude is underestimated when compared to observations found at the leaf scale. This finding raises the question of whether the current parameterization of the stomatal control is sufficient to correctly reproducing the 385 leaf's behavior in nature. A first approach in this regard could be the implementation of a process that artificially regulates the stomatal resistance such that the plant can stay close to its temperature optimum. Furthermore, these temperature optima can adapt to different climatological conditions (Berry and Bjorkman, 1980), which is a crucial aspect in a future warmer climate.
In this sense, it would be valuable to further study the phenomenon of leaf thermoregulation, because in a warming climate many plants may come close to their temperature limits. Thus, they would increasingly be dependent on the cooling effect due 390 to the leaf's ability to regulate its temperature. Therefore, estimates of carbon uptake by the biosphere may crucially depend on the accurate representation of leaf thermoregulation. technique (Ershadi et al., 2014). Overall, the simulated fluxes of Classic, as well as of CEBa, fit the observations quite well.

415
At the tropical site, Classic underestimates the daytime maximum in radiative temperature, while CEBa overestimates it for Tharandt. In general, the Classic scheme shows numerical induced fluctuations in the heat fluxes even in the average over multiple years, while CEBa exhibits smooth curves such as seen in the observations. For the tropics, the Classic scheme shows in this regard an unusual phase shift in the sensible heat flux compared to the observations and the CEBa scheme. Since the radiative temperature is not phase-shifted, it could be caused by the abovementioned numerical instabilities rather than by 420 incorrectly simulated heat storages. For Tharandt, both schemes show a significant difference in the Bowen ratio: the latent heat flux is overestimated, and the sensible heat flux underestimated. The cause of this bias could not be explained in this study and should be further investigated in the future.
In addition to the offline single-site experiments, we evaluated the performance of CEBa with respect to the simulated 2m-temperature in a global coupled AMIP experiment over 30 years (from 1979-2008). We find a similar pattern in the 425 difference of the simulated near-surface temperature between CEBa and Classic as in the difference between SkIn + and Classic (see Fig. 5 in Heidkamp et al., 2018). Only a yet unexplained slight temperature increase in the high latitudes (mainly over Siberia) remains. That also leads to a significant bias reduction in the near-surface temperatures simulated by the CEBa scheme.
However, despite the potentially more realistic representation of the processes in the canopy layer, the CEBa scheme does not achieve the same model bias reduction (7.4 %) as the SkIn + scheme (8.8 %). 430 We find that the new CEBa scheme correctly reproduces the diurnal cycle of T sfc , H, and LE at the site-level for a temperate as well as for a tropical forest. Only for the daytime maximum in radiative temperature in Tharandt, CEBa shows a larger temperature bias than JSBACH Classic, but it removes the unrealistic drop in sensible heat storage during the afternoon at the tropical forest site. Regarding the performance of the CEBa scheme on regional and global scales, we conclude that CEBa leads to a significant bias reduction in simulated near-surface temperatures compared to JSBACH Classic, but has a slightly 435 weaker performance than SkIn + . To improve the results of the CEBa scheme, an extensive tuning and technical effort would be necessary, which is beyond the scope of this study. However, in general, the CEBa scheme provides, due to the more realistic treatment of energy exchange and its larger number of degrees of freedom, the possibility of a deeper understanding of the complex processes in the canopy layer.

440
For the estimation of u cas , we use an approach according to Xue et al. (1991). The basic idea is the concept of a so-called transition height z trans . This height can be either below or above the reference height, depending on the height of the roughness layer, which is determined by the height of the vegetation. Figure B1 illustrates the concept considering two different cases for the transition height. For tall vegetation, the first case, the transition height exceeds the reference height of the LAL z ref .
For shallow vegetation, the transition height may be lower than the LAL. In this second case, the increase in wind speed of the 445 logarithmic wind profile is scaled by an adjustment factor (G 2 = 0.75). Therefore, the wind speed in the CAS at the top of the vegetation z veg can be calculated as

Roughness layer
Prandtl layer 2 nd Case z trans < z ref

Appendix C: Physics of leaf thermoregulation in CEBa
In the context of the CEBa scheme, the negative correlation between the leaf temperature excess and the air temperature can be explained by the energy balance of the leaf as follows: Assuming low temperatures and constant net radiation, an increase in air temperature will initially raise the leaf temperature to maintain the energy equilibrium by keeping a constant temperature 455 difference and sensible heat flux. However, while the leaf temperature increases linearly, the saturated specific humidity of the leaf q sat (T leaf ) increases exponentially following Clausius-Clapeyron theory. Thus, at higher temperatures the latent heat flux increasingly outweighs the sensible heat flux due to transpiration. Above a certain air temperature, the latent heat flux of the Figure C1. Example case for the analytically solved leaf energy balance equation: Negative correlation between leaf temperature excess and air temperature as well as three different cases of the energy balance (Rnet in violet, H leaf in red and LE leaf in blue) according to Eq.

Simple Energy Balance Model of a Leaf
(C4).
leaf exceeds the available energy flux, resulting in a sensible heat flux towards the leaf (Fig. C1). This behavior is known as the oasis effect and has also been found for vegetated canopies (Taha et al., 1991).

460
To explain this process quantitatively, we consider the energy balance of the leaf and solve it analytically for the temperature difference T leaf − T air . Ignoring heat storages within the leaf, the energy balance can be written by where R net is the net radiation, LE leaf the latent-and H leaf the sensible heat flux of the leaf. The latent and sensible heat flux can be expressed as where q air is the specific humidity of the air, r b is the aerodynamic leaf resistance and r s the stomatal resistance. The terms c 1 = ρL v as well as c 2 = ρc p can be seen as constant, in first order. Here ρ denotes the air density, L v the enthalpy of vaporization and c p the specific heat capacity of air at constant pressure. The saturated specific humidity of the latent heat fluxes can be linearized at T = T air as 470 LE ≈ c 1 q sat (T air ) − q air r b + r s + c 1 α T leaf − T air r b + r s where α is the slope of the saturation specific humidity at T air . With the abovementioned approximations, the leaf energy balance can be solved analytically and written in terms of the leaf temperature excess, as follows with c 3 = 1 + αc 1 /c 2 . With this equation, we can quantify the atmospheric conditions and leaf properties that lead to an 475 equilibrium between the leaf and the air temperature (equilibrium point). Equation C4 shows that an increase in net radiation or relative humidity decreases the leaf temperature excess and the plant's ability to buffer its temperature against air temperature.
Author contributions. All authors designed the experiments. MH performed and analyzed the simulations, developed the model code, and prepared the paper. All authors contributed to generating ideas, writing the paper, and discussing the results.
Competing interests. The authors declare that they have no conflict of interest.

480
Acknowledgements. The study was supported by the Max Planck Society for the advancement of science. The use of the supercomputer facilities at the Deutsches Klimarechenzentrum (DKRZ) is acknowledged. We would like to thank all people contributing to the FLUXNET network in particular these involved in the creation of the FLUXNET2015 Dataset (Pastorello et al., 2020). Especially, we thank Christian Bernhofer, Thomas Grünwald, Uta Moderow, Markus Hehn, Uwe Eichelmann, and Heiko Prasse for compiling the Tharandt dataset. The work of the tower team of the tropical FLUXNET site around Mike Goulden is also gratefully acknowledged.

485
The article processing charges for this open-access publication were covered by the Max Planck Society.