# Increased control of vegetation on global terrestrial energy fluxes

Giovanni Forzieri, Diego Miralles, Philippe Ciais, Ramdane Alkama, Youngryel Ryu, Gregory Duveiller, Ke Zhang, Eddy Robertson, Markus Kautz, Brecht Martens, et al.

► To cite this version:

Giovanni Forzieri, Diego Miralles, Philippe Ciais, Ramdane Alkama, Youngryel Ryu, et al. Increased control of vegetation on global terrestrial energy fluxes. Nature Climate Change, 2020, 10 (4), pp.356-362. 10.1038/s41558-020-0717-0. hal-02968000

HAL Id: hal-02968000

https://hal.science/hal-02968000v1

Submitted on 2 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Increased control of vegetation on global
terrestrial energy fluxes

Giovanni Forzieri ID <fnref n="1" />✉, Diego G. Miralles ID ², Philippe Ciais ID ³, Ramdane Alkama ID <fnref n="1" />, Youngryel Ryu ⁴,
Gregory Duveiller ID <fnref n="1" />, Ke Zhang ID ⁵, Eddy Robertson ⁶, Markus Kautz ID ⁷, Brecht Martens ²,
Chongya Jiang ID ⁴, Almut Arneth ⁸, Goran Georgievski ⁹,¹⁰, Wei Li ID ¹¹, Guido Ceccherini ID <fnref n="1" />,
Peter Anthoni ⁸, Peter Lawrence ID ¹², Andy Wiltshire ⁶, Julia Pongratz ID <fnref n="13" />,¹⁴, Shilong Piao ID ¹⁵,
Stephen Sitch ID ¹⁶, Daniel S. Goll ID ¹⁷, Vivek K. Arora ¹⁸, Sebastian Lienert ID ¹⁹, Danica Lombardozzi ID ¹²,
Etsushi Kato ID ²⁰, Julia E. M. S. Nabel ID ¹⁴, Hanqin Tian ID ²¹, Pierre Friedlingstein ID <fnref n="22" /> and
Alessandro Cescatti<fnref n="1" />

Changes in vegetation structure are expected to influence the redistribution of heat and moisture; however, how variations in the leaf area index (LAI) affect this global energy partitioning is not yet quantified. Here, we estimate that a unit change in LAI leads to 3.66 ± 0.45 and -3.26 ± 0.41 W m⁻² in latent (LE) and sensible (H) fluxes, respectively, over the 1982-2016 period. Analysis of an ensemble of data-driven products shows that these sensitivities increase by about 20% over the observational period, prominently in regions with a limited water supply, probably because of an increased transpiration/evaporation ratio. Global greening has caused a decrease in the Bowen ratio (B = H/LE) of -0.010 ± 0.002 per decade, which is attributable to the increased evaporative surface. Such a direct LAI effect on energy fluxes is largely modulated by plant functional types (PFTs) and background climate conditions. Land surface models (LSMs) misrepresent this vegetation control, possibly due to underestimation of the biophysical responses to changes in the water availability and poor representation of LAI dynamics.

The surface energy partitioning, which results from the con-
version of available energy into LE and H, exerts a strong
control on the state of the atmospheric boundary layer, the
lowest layer of the troposphere that is in contact with the surface
of the Earth. It propagates changes in land surface properties to the
atmosphere<fnref n="1" />,² by regulating land-atmosphere feedback effects and
influencing the global cycles of water and energy³. Whereas most
research has been devoted to exploring the effects of soil moisture
on surface energy partitioning⁴,⁵, vegetation density can also have
an important role in the modulation of the surface energy budget⁶⁻⁸.
In fact, changes in vegetation structure and physiology associated
with the ongoing global and persistent increase in LAI⁹ are expected
to influence canopy conductance, aerodynamic properties and the
albedo of ecosystems, ultimately affecting water and energy fluxes
between land and atmosphere¹⁰,¹¹. Consequent variations in cli-
mate modulate the interplay between LAI-related biophysical pro-
cesses and the surface energy partitioning¹²,<fnref n="13" />. Rising atmospheric
CO₂ concentrations further affect the vegetation control on sur-
face energy partitioning by reducing stomatal conductance¹⁴ and

therefore transpiration per unit of leaf area, ultimately leading to an
increasing ratio of carbon gain to water loss (water-use efficiency)¹⁵.
Although the influence of global greening feedback effects on sur-
face temperature has been recognized and assessed¹²,¹⁶, its more
direct impacts on the surface energy partitioning at the planetary
level have not been explicitly explored.

Disentangling the role of LAI from the contribution of other direct drivers on the surface energy partitioning is challenging due to the variety of land-atmosphere interactions that occur across multiple spatial and temporal scales¹⁷. For this scope, the availability of observations from field experiments¹⁸ and flux tower data¹⁹ is limited in terms of the number of stations and geographical coverage, enabling only a partial characterization of the spatiotemporal variability of the phenomena. On the modelling side, LSMs—the land component of Earth system models used to predict future climate trajectories—include LAI as a key prognostic variable and its interactions with surface biophysics, hydrology and biogeochemistry are represented through equations of varying complexity²⁰. However, LSMs show important limitations in reproducing the interplay

<footnote n="1">¹European Commission, Joint Research Centre, Ispra, Italy. ²Hydro-Climate Extremes Lab (H-CEL), Department of Environment, Ghent University, Ghent, Belgium. ³LSCE CEA CNRS UPSACLAY IPSL, Gif-sur-Yvette, France. ⁴Department of Landscape Architecture and Rural Systems Engineering, Seoul National University, Seoul, South Korea. ⁵College of Hydrology and Water Resources, Hohai University, Nanjing, China. ⁶Met Office, Exeter, UK. ⁷Forest Research Institute Baden-Württemberg, Freiburg, Germany. ⁸KIT, IMK-IFU, Garmisch-Partenkirchen, Germany. ⁹Helmholtz-Zentrum Geesthacht, Institute of Coastal Research, Geesthacht, Germany. ¹⁰Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Bremerhaven, Germany. ¹¹Ministry of Education Key Laboratory for Earth System Modeling, Department of Earth System Science, Tsinghua University, Beijing, China. ¹²NCAR, Boulder, CO, USA.</footnote>

<footnote n="13">¹³Ludwig-Maximilians-Universität München, Munich, Germany. ¹⁴Max Planck Institute for Meteorology, Hamburg, Germany. ¹⁵Department of Ecology, Peking University, Beijing, China. ¹⁶College of Life and Environmental Sciences, University of Exeter, Exeter, UK. ¹⁷Department of Geography, University of Augsburg, Augsburg, Germany. ¹⁸Canadian Centre for Climate Modelling and Analysis, Environment Canada, University of Victoria, Victoria, British Columbia, Canada. ¹⁹Climate and Environmental Physics, Physics Institute and Oeschger Centre for Climate Change Research, University of Bern, Bern, Switzerland. ²⁰Institute of Applied Energy (IAE), Tokyo, Japan. ²¹School of Forestry and Wildlife Sciences, Auburn University, Auburn, AL, USA.</footnote>

<footnote n="22">²²College of Engineering, Mathematics and Physical Sciences, University of Exeter, Exeter, UK. ²³e-mail: giovanni.forzieri@ec.europa.eu</footnote>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 1 |** Sensitivity of surface energy partitioning to LAI changes. **a,** Sensitivity of latent heat ($\frac{\partial LE}{\partial LAI}$) to changes in LAI computed for the 1982–2016 period and binned as a function of climatological mean *P* and *T*. Black dots indicate bins with average values that are statistically different from zero (Student's *t*-test; *P*-value ≤ 0.05). **b,** Temporal variations in ($\frac{\partial LE}{\partial LAI}$) (ensemble average ± s.e.) computed over a 13-year moving window for moisture supply-limited and atmospheric demand-limited regions, and the whole globe. **c,** Temporal variations in ($\frac{\partial LE}{\partial LAI}$) computed separately over the 1982–1999 and 2000–2016 periods for different regions and expressed in terms of absolute variations ($\Delta_{abs}$) and relative variations ($\Delta_{rel}$). Results for each single LAI-ET combination are shown using different symbols; those with a black outline represent ensemble averages (both computed for each ET product and for the whole set of combinations labelled as 'Ensemble' in the key); black dots indicate statistically significant changes in sensitivity (Student's *t*-test, *P*-value ≤ 0.05). The spatial domains of supply- and demand-limited regions are shown in Supplementary Fig. 5. **d–i,** Data as in **a–c;** the sensitivity of sensible heat ($\frac{\partial H}{\partial LAI}$) and the Bowen ratio ($\frac{\partial B}{\partial LAI}$) to changes in LAI are shown in **d–f and g–i,** respectively.</figcaption></figure>

between vegetation and climate due to an incomplete understanding and model representation of biophysical processes<sup>21,22</sup>. These drawbacks inevitably hamper our understanding of land-atmosphere interactions based on in situ observations and model predictions. Conversely, the increasing availability and accessibility of satellite remote-sensing products that address the physical state of the land surface may overcome these limitations and offer robust global datasets for model evaluation and process understanding.

Here, we investigate the impacts of greening on LE, *H* and the Bowen ratio for 1982–2016 at the global scale, using four observation-driven products of evapotranspiration (ET) and climate drivers in combination with three long-term satellite LAI datasets. We refer to the growing-season-averaged LAI as a diagnostic variable of vegetation density. The sensitivity of LE, *H* and *B* to changes in LAI is quantified as partial derivatives from multiple linear regressions ($\partial Z/\partial LAI$ where *Z* is any of the energy terms). The potentially confounding direct effect (not through LAI) of climate drivers such

as precipitation and temperature on LE, *H* and *B* has been factored out by considering these variables among the predictors of the linear model (Methods). Furthermore, interannual variations of both response variable and predictors have been used in the regression to rule out possible long-term dependencies between covariates. Long-term effects in surface energy terms that are attributable strictly to greening ($\delta Z^{L-ALI}$) are then quantified by combining sensitivity estimates with long-term trends in LAI. Sensitivity ($\partial Z/\partial LAI$) and effects ($\delta Z^{L-ALI}$) are derived for each combination of LAI and energy-flux datasets resulting in a 12-member ensemble of observation-based estimates. To account for differences across datasets<sup>23</sup>, the ensemble average is calculated and the corresponding standard error (s.e.) is obtained (Methods). Finally, we compared the metrics derived from satellite-based observations with those computed on factorial simulations of ten state-of-the-art LSMs, to assess the ability of models to represent the interplay between vegetation changes and surface energy partitioning (Methods).
a

<figure><img src="image_12.png" /></figure>

b

<figure><img src="image_13.png" /></figure>

c

<figure><img src="image_14.png" /></figure>

d

<figure><img src="image_15.png" /></figure>

e

<figure><img src="image_16.png" /></figure>

f

<figure><img src="image_17.png" /></figure>

g

<figure><img src="image_18.png" /></figure>

h

<figure><img src="image_19.png" /></figure>

i

<figure><img src="image_20.png" /></figure>

**Fig. 2** | Changes in surface energy partitioning associated with long-term trends in LAI. **a–c**, Spatial pattern (**a**), seasonal variability (**b**) and climate space (**c**) of LAI-related trend in latent heat (δLE<sup>LAI</sup>; W m<sup>-2</sup> per decade) computed for the 1982–2016 period. **d–i**, Data as in **a–c**; the LAI-related trends in sensible heat (δH<sup>LAI</sup>; W m<sup>-2</sup> per decade) and Bowen ratio (δB<sup>LAI</sup>; per decade) are shown in **d–f** and **g–i**, respectively. Black dots in **a,d,g** show pixels for which both the ensemble average LAI trend and sensitivity are significant (Mann–Kendall test and Student's t-test, respectively; *P*-value < 0.05). Values in **b,e,h** are binned as a function of time and latitudinal gradient, and black dots show bins with average values that are statistically different from zero (Student's t-test; *P*-value < 0.05). Values in **c,f,i** are binned as a function of climatological mean *P* and *T* and black dots are as in **b,e,h**.

Increased sensitivity of surface energy partitioning on LAI

Estimates of the sensitivity of the energy partitioning terms to LAI (equation (1), Methods) quantified for the 1982–2016 period show a clear dependence on the background climate (Fig. 1a,d,g). In accordance with previous model-based studies⁶, the increase in LAI enhances LE globally (3.66 ± 0.45 W m⁻² per unit of leaf area; Supplementary Table 1) and particularly in warm and dry regions, as a consequence of the increase in evaporative surface (Fig. 1a). In these regions, despite potential limitations of soil moisture on LE due to low rainfall, at the interannual timescale an increase in LAI is associated with an increase in LE through complex adjustments of the LAI sustaining LE, such as root development, access to groundwater²⁴ and phenological seasonal shifts²⁵. Given that LE and *H* represent competitive pathways for energy release from the land surface, *H* shows opposite patterns of sensitivity to LAI than LE, with an average negative sensitivity of −3.26 ± 0.41 W m⁻² per unit of leaf area (Fig. 1d, Supplementary Table 1). Ultimately, changes in *B* are inversely related to LAI (−0.14 ± 0.02 per unit of leaf area; Fig. 1g and Supplementary Table 1), as the increase in leaf area favours the dissipation of available energy by evaporating water, leading to surface cooling and a subsequent reduction in *H*¹²,¹⁶. These findings

emphasize the importance of interannual vegetation controls on climate, particularly during extreme events such as meteorological droughts and heatwaves, when a higher LAI can effectively dampen the increases in land surface temperature by evaporative cooling, although at the expense of further drying out the soils²⁶.

Exploring the temporal variation in the sensitivities using moving windows with different lengths (for brevity, a 13-year window is shown here), we found substantial changes over the observational period, particularly from 2000 onwards (Fig. 1b,e,h). Globally, we found a significant (*P*-value ≤ 0.05) relative increase of 20–24% (Δrel) in the value of the sensitivity of energy fluxes to LAI over 2000–2016 compared with 1982–1999 (grey circles in Fig. 1c,f,i), suggesting an increasing control of energy fluxes from terrestrial vegetation. Despite the relevant spread observed across single LAI products, global trends are largely consistent (Extended Data Fig. 1) and show a dependence of the variation in sensitivity to the changes in interannual LAI variability (Supplementary Fig. 1). Although we cannot exclude a possible contamination of the variability in LAI from the temporal variations in satellite platforms and sampling density, we stress that climate variations are very likely to have a major role on the emerging signal. This finding is supported by
<figure><img src="image_21.png" /><figcaption>**Fig. 3** | Comparison of LAI and climate effects on surface energy partitioning. **a–c**, Effects (x) of the long-term trends in LAI, T, P and SW<sub>in</sub> on the latent heat (δLE<sup>+</sup>) (**a**), sensible heat (δH<sup>+</sup>) (**b**) and Bowen ratio (δB<sup>+</sup>) (**c**) for moisture supply-limited and atmospheric demand-limited regions and the whole globe. Box plots represent the 12-member ensemble of observation-driven products. The spatial domains of supply-limited and demand-limited regions are shown in Supplementary Fig. 5. **d–f**, Density plot of pixel values of overall trends in latent heat (δLE) (**d**), sensible heat (δH) (**e**) and Bowen ratio (δB) (**f**) versus the corresponding LAI-related effect δLE<sup>LAI</sup>, δH<sup>LAI</sup> and δB<sup>LAI</sup>, respectively. Red labels indicate the fraction of the global domain that falls within each quadrant delineated by red lines; blue circles show the global-scale estimates with labels referring to the fraction of the overall trend in the energy term explained by the effect of LAI alone.</figcaption></figure>

the analysis of the gradients in sensitivity across space (Extended Data Fig. 2, Methods). The dependence of the sensitivity on the aridity index over the two observation periods is statistically identical, despite the change in the observation system (from AVHRR to MODIS). These results further corroborate the relevance of climate change on the temporal variation in the sensitivity.

The possible mechanisms responsible for such emerging variations in the sensitivity were investigated by disaggregating the signal for regions in which evaporation is limited by atmospheric demand or by water supply (Methods). We found that the largest absolute variations occur in regions limited by the supply of moisture in which the change in the sensitivity is fourfold higher than in demand-limited regions (for example, the Δ<sub>abs</sub> of sensitivity of LE is around 1.8 and 0.4 W m<sup>-2</sup>, respectively; Fig. 1c). The concomitant strong positive trends in temperature and moderately negative trends in precipitation lead to a progressive transition to warmer and drier conditions (Extended Data Fig. 3c–f), which is represented by a shift towards the top-left corner in Fig. 1a,d,g. In moisture supply-limited environments, such changes in climate background have probably increased the ratio between transpiration and evaporation, leading to an enhanced biological control on evapotranspiration. This hypothesis is supported by the widespread increasing trend in the fraction of transpired water to the total evapotranspiration (Extended Data Fig. 3i,j) that has also been documented in previous studies<sup>27</sup>. The recent climate-induced expansion of areas limited by water supply has presumably amplified this process (Extended Data Fig. 4a).

We found that, consistently with expectations, temporal changes in sensitivity are lower for the datasets that explicitly account for the direct CO₂ effects on stomatal conductance and transpiration (PLSH and BESS) than for the other products (GLEAM and MTE). However, a significant increase in sensitivity (Δ<sub>rel</sub> > 10%) emerges even for ET products that consider CO₂ effects (Fig. 1c,f,i and Extended Data Fig. 1), therefore suggesting that the increase in water-use efficiency<sup>28</sup> cannot fully offset the emerging climate signal of increasing control of plant leaves on terrestrial energy fluxes.

## Greening has a key role in surface energy partitioning

The effects of greening on surface energy partitioning are derived by multiplying the observed sensitivity by the long-term trend in growing-season-averaged LAI as quantified for the 1982–2016 period (Extended Data Figs. 3a,b and 5a–c)—that is, by applying the previously described methodology<sup>12</sup> (equation (4), Methods). Results show that the variations in LAI that have occurred over the past three-and-a-half decades as such have led to a significant increase in LE across a large part of the globe (Fig. 2a,c), particularly in moisture supply-limited regions (0.41 ± 0.09 W m<sup>-2</sup> per decade; Supplementary Table 1). Such a pronounced effect of greening results from a combination of moderately positive trends in LAI (Extended Data Fig. 3a,b) and the high sensitivity of the latent heat fluxes to LAI in those regions (Fig. 1a), consistently with previous findings<sup>27,29</sup>. By contrast, atmospheric demand-limited regions show a limited effect of changes in LAI on LE trends (Fig. 2a,c), primarily due to the low sensitivity of evapotranspiration to changes
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

**Fig. 4** | Contribution of different PFTs to the LAI control on energy partitioning. **a**, Spatial map of cover fractions of PFTs (forests, natural grasses and croplands). **b**, Relative contribution of each PFT to the global land area (Area), long-term variations in LAI (δLAI), sensitivity of energy partitioning terms to changes in LAI (δZ/δLAI), and long-term effects in surface energy terms due to long-term variations in LAI (δZ<sup>LAI</sup>), where Z is LE, H or B. **c–e**, Cover fractions of PFTs binned as a function of the long-term variations in LAI (δLAI) and the sensitivity of latent heat (δLE/δLAI) (**c**), sensible heat (δH/δLAI) (**d**) and Bowen ratio (δB/δLAI) (**e**) to LAI changes. Average values are shown for each PFT in circles.

in LAI in these areas (Fig. 1a). As expected, trends in *H* associated with greening are opposite to those of LE, due to their reversed sensitivity (Fig. 2d,f). The combination of increasing LE and decreasing *H* trends attributable to the greening signal led to a widespread decrease in *B* (−0.010 ± 0.002 per decade; Fig. 2g,i, Supplementary Table 1). As these effects of LAI are larger in water-limited regions associated with high *B* values, the recent greening may have reduced the spatial variability of surface energy partitioning across Earth, ultimately affecting the strength of land–atmosphere coupling and the dynamics of the boundary layer.

The seasonality of the sensitivity shows peak values at the onset of the growing season, particularly in temperate regions in the Northern Hemisphere (Fig. 2b,e,h). This suggests that, during the growing season, the partition of surface-available energy is more closely controlled by variations in the LAI, mainly due to an increase in the ratio of transpiration to evaporation. Climate-driven changes in plant phenology may further amplify these mechanisms, particularly at the beginning of the growing season due to the expected increase in transpiration associated with an earlier onset of vegetation green-up³⁰. Furthermore, we note that the LAI-related trend in energy partitioning shows a modest seasonal pattern also at high latitudes with a change in sign that is more evident for LE (Fig. 2b,e). In this climate zone, the interplay between LAI and energy partitioning seems to be dominated by radiative terms during cold seasons (a LAI-related reduction in albedo increases *H* and decreases LE during conditions of snow cover), but is still dominated by changes in the evaporative surface during the warm seasons (strong positive control of LAI on transpiration), consistent with previous studies³¹.

The methodology used to quantify the effects of the long-term trends in LAI on surface energy partitioning at the annual level is then applied to each predictor of the regression to assess their relative contributions. A large variability across supply- and demand-limited zones emerges for the effects associated with LAI, air temperature (T), precipitation (P) and short-wave incoming radiation (SW<sub>in</sub>) (Fig. 3a–c), reflecting the spatial variations in long-term trends and sensitivity of each variable (Extended Data Figs. 3 and 5). Although our assessment did not include possible interactions among drivers, it unequivocally shows that LAI has a larger amount of control than the direct effects of T, P and radiation on the trends in energy fluxes (Fig. 3a–c). Notably, the low contribution of P mostly results from its low and spatially variable long-term trend (Extended Data Fig. 3e,f). In addition to the comparison among single drivers, we found that the effects of LAI are concordant in sign with the overall trends in the energy partitioning terms (Supplementary Fig. 2) for more than 63% of the vegetated land (Fig. 3d–f) and explain a considerable fraction of their variance (45–63%; Fig. 3d–f). These findings emphasize the importance of LAI trends in affecting the long-term variations in surface energy partitioning, in particular by amplifying the release of energy as latent heat (Fig. 2a,d,g and Supplementary Fig. 2).

## Leaf control and PFTs

The spatial distribution of PFTs (Fig. 4a) modulates the effects of LAI changes on energy partitioning. In the observation period, forecasts show a strong increase in LAI, which is predominantly driven by climate change and CO₂ fertilization⁹, and provide the largest
<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

**Fig. 5** | Comparison of observational and LSM results. **a**, Temporal variations in the sensitivity of LE to LAI changes ($\frac{\partial LE}{\partial LAI}$) retrieved from an ensemble of 10 LSMs (TRENDY v.7) under the S3 scenario (changes in CO₂, climate and land use) and computed over 13-year moving windows for supply- and demand-limited regions and the whole globe (ensemble average ± s.e.). Labels at the bottom of the panel report the relative changes in sensitivities between the 1982-1999 period and the 2000-2016 period ($\Delta_{e}$). Significance was calculated by Student's t-test; \* *P*-value ≤ 0.05. The spatial domains of supply- and demand-limited regions are shown in Supplementary Fig. 5. **b**, As in **a**, but for the S1 scenario (changes in CO₂ only). **c**, Spatial patterns of the differences in ($\frac{\partial LE}{\partial LAI}$) between LSMs and satellite-driven products (TRENDY S3 minus SAT) computed for the 1982-2016 period as the median of all data-model combinations. **d**, Values in **c** are binned as a function of the climatological mean of P and T. **e**, Values in **c** are binned as a function of satellite and modelled absolute interannual variations in LAI. Black dots in **c**–**e** show pixels and bins with average values that are statistically different from zero (Student's t-test; *P*-value < 0.05).

contribution to the global signal of greening (48%; Fig. 4b). However, they are typically characterized by a low sensitivity of the energy terms to LAI that is possibly due to a more-conservative and even use of water resources supported by a deeper rooting system³² and by their abundance in demand-limited regions. Therefore, when the two terms are combined (greening and sensitivity; equation (4), Methods), forests contribute 21–27% to the global effect ($\delta Z^{L}_{AI}$, where Z is any of the energy terms). By contrast, natural grasses show lower greening rates compared with forests but have a larger sensitivity driven by the rapid dynamics of LAI and shallow soil moisture, which are typical of these ecosystems. This combination results in a contribution of 32–38% to the global effect of greening on the surface energy partitioning. Croplands, despite their limited fractional cover (25%), have a more-important role in affecting the surface energy partitioning by contributing 41–43% to the global signal. This derives from the combination of relatively high values of both sensitivity and greening (Fig. 4c–e), the latter of which is possibly driven by agricultural intensification, which occurred in many regions of the world during the past decades³³.

## Possible sensitivity bias in LSMs

The ability of LSMs to reproduce the interplay between LAI and surface energy partitioning has important implications for the prediction of future land-climate interactions²². To assess this ability, we replicated the analysis performed on observation-driven products on an ensemble of 10 state-of-the-art LSMs (TRENDY v.7)³⁴. Model simulations show that the CO₂ effect on stomatal conductance substantially reduces the sensitivity of LE to changes in LAI (–11%; Fig. 5b, S1 scenario: changes in CO₂ only) to the point of offsetting the modest but significant signal that originated from changes in climate (8%; Extended Data Fig. 6c11, S3–S1 scenario: changes in

climate and land use only). This ultimately leads to a net signal of no change in sensitivity in the scenario in which all factors are varied (Fig. 5a, S3 scenario: changes in CO₂, climate and land use). This pattern, in combination with the predicted reduction in supply-limited zones (Extended Data Fig. 4b), leads to an overall decrease in the coupling between LE and water availability. These simulated patterns are in clear contrast to those retrieved from observation-driven products (Fig. 1 and Extended Data Fig. 4a). Such divergence suggests that LSMs may overestimate the sensitivity to CO₂ and underestimate the biophysical response of ecosystems to changes in water availability. This pattern emerges also from the systematic model underestimation of the fraction of transpired water to the total evapotranspiration³⁵. Under a scenario of warming, this bias of LSMs could ultimately lead to an underestimation of summer droughts sustained by anticipated spring phenology³⁶. In addition, the large spread across LSMs (Fig. 5a,b and Extended Data Fig. 6) highlights the large structural uncertainty in the model representation of the phenomena. In fact, even if the ensemble is driven by a common climate forcing, the model structure and parameterization show a large effect on the energy partitioning⁴⁷.

Overall, focusing on the sensitivity derived over the whole period and assuming the observation-based sensitivity of LE to LAI as reference, the tested LSMs show an overestimation of the sensitivity for tropical and boreal zones, and an underestimation over arid-temperate zones (average data-model discrepancy of 0.89±0.44 and 2.9±0.53 W m⁻² per unit LAI change, respectively) (Fig. 5c,d). Such differences are associated with the underestimation and overestimation of ΔLAI in the two regions, respectively (Fig. 5e). Models show sensitivities and LAI trends with a bias of an opposite sign with respect to observational retrievals for about half of the globe and such compensatory effects hide the effective
data-model discrepancies in the resulting global effect of greening on the partitioning of the surface energy fluxes (Extended Data Fig. 7). Although disparities with respect to observation-based findings are conditioned by the accuracy of satellite retrievals—which are particularly important in the tropics where LAI estimations tends to saturate—these results emphasize the current uncertainty generated by the approximate model representation of key vegetation-mediated biophysical processes.

## Conclusions

Our analysis provides observational evidence that changes in vegetation density (LAI) during the past three-and-a-half decades have had an important role in the surface energy partitioning, by favouring the release of energy through LE rather than H. This increased control of vegetation appears to be plausibly connected to the exacerbation of water-limited conditions and the progressive increase in evaporative surface associated with the global greening, and occurs despite the counteracting effect of CO₂ fertilization on stomatal conductance. As land feedback effects on climate are linked to vegetation status and activity, future land geoengineering could play an important part in modulating the strength of surface forcing⁹. Furthermore, our results reinforce the importance of considering the co-variability of soil moisture and vegetation dynamics for the effective appraisal of land-atmosphere coupling (for which usually the focus was exclusively on soil moisture variability and patterns)⁷, particularly in view of the expected increase in LAI³⁸ and drought conditions³⁹ across most of the globe. Finally, our data-model comparison emphasizes the need to better account for the influence of vegetation changes in energy partitioning to improve climate model projections. Fostering model representation of vegetation-atmosphere interactions with observation-driven estimates will ultimately enhance the reliability of future climate predictions.

## Online content

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41558-020-0717-0.

Received: 17 September 2018; Accepted: 3 February 2020;
Published online: 16 March 2020

## References

1. Findell, K. L., Gentine, P., Lintner, B. R. & Kerr, C. Probability of afternoon precipitation in eastern United States and Mexico enhanced by high evaporation. *Nat. Geosci.* **4**, 434–439 (2011).

2. Pielke, R. A. Influence of the spatial distribution of vegetation and soils on the prediction of cumulus convective rainfall. *Rev. Geophys.* **39**, 151–177 (2001).

3. Pitman, A. J. The evolution of, and revolution in, land surface schemes designed for climate models. *Int. J. Climatol.* **23**, 479–510 (2003).

4. Seneviratne, S. I. et al. Investigating soil moisture-climate interactions in a changing climate: a review. *Earth Sci. Rev.* **99**, 125–161 (2010).

5. Dirmeyer, P. A. The terrestrial segment of soil moisture-climate coupling. *Geophys. Res. Lett.* **38**, L16702 (2011).

6. Puma, M. J., Koster, R. D. & Cook, B. I. Phenological versus meteorological controls on land-atmosphere water and carbon fluxes. *J. Geophys. Res. Biogeosci.* **118**, 14–29 (2013).

7. Williams, I. N. & Torn, M. S. Vegetation controls on surface heat flux partitioning, and land-atmosphere coupling. *Geophys. Res. Lett.* **42**, 9416–9424 (2015).

8. Bagley, J. E. et al. The influence of land cover on surface energy partitioning and evaporative fraction regimes in the U.S. Southern Great Plains. *J. Geophys. Res. Atmos.* **122**, 5793–5807 (2017).

9. Zhu, Z. et al. Greening of the Earth and its drivers. *Nat. Clim. Change* **6**, 791–795 (2016).

10. Anderson, R. G. et al. Biophysical considerations in forestry for climate protection. *Front. Ecol. Environ.* **9**, 174–182 (2011).

11. Bonan, G. B. Forests and climate change: forcings, feedbacks, and the climate benefits of forests. *Science* **320**, 1444–1449 (2008).

12. Forzieri, G., Alkama, R., Miralles, D. G. & Cescatti, A. Satellites reveal contrasting responses of regional climate to the widespread greening of Earth. *Science* **356**, 1180–1184 (2017).

13. Pitman, A. J. et al. Importance of background climate in determining impact of land-cover change on regional climate. *Nat. Clim. Change* **1**, 472–475 (2011).

14. Berry, J. A., Beerling, D. J. & Franks, P. J. Stomata: key players in the earth system, past and present. *Curr. Opin. Plant Biol.* **13**, 232–239 (2010).

15. Keenan, T. F. et al. Increase in forest water-use efficiency as atmospheric carbon dioxide concentrations rise. *Nature* **499**, 324–327 (2013).

16. Zeng, Z. et al. Climate mitigation from vegetation biophysical feedbacks during the past three decades. *Nat. Clim. Change* **7**, 432–436 (2017).

17. Wilson, K. B. et al. Energy partitioning between latent and sensible heat flux during the warm season at FLUXNET sites. *Water Resour. Res.* **38**, 1294 (2002).

18. Sellers, P. J. et al. BOREAS in 1997: experiment overview, scientific results, and future directions. *J. Geophys. Res. Atmos.* **102**, 28731–28769 (1997).

19. Baldocchi, D. et al. FLUXNET: a new tool to study the temporal and spatial variability of ecosystem-scale carbon dioxide, water vapor, and energy flux densities. *Bull. Am. Meteorol. Soc.* **82**, 2415–2434 (2001).

20. Piao, S. et al. Evaluation of terrestrial carbon cycle models for their response to climate variability and to CO₂ trends. *Glob. Change Biol.* **19**, 2117–2132 (2013).

21. Duveiller, G. et al. Biophysics and vegetation cover change: a process-based evaluation framework for confronting land surface models with satellite observations. *Earth Syst. Sci. Data* **10**, 1265–1279 (2018).

22. Forzieri, G. et al. Evaluating the interplay between biophysical processes and leaf area changes in land surface models. *J. Adv. Model. Earth Syst.* **10**, 1102–1126 (2018).

23. Jiang, C. et al. Inconsistencies of interannual variability and trends in long-term satellite leaf area index products. *Glob. Change Biol.* **23**, 4133–4146 (2017).

24. Fan, Y., Li, H. & Miguez-Macho, G. Global patterns of groundwater table depth. *Science* **339**, 940–943 (2013).

25. Morisette, J. T. et al. Tracking the rhythm of the seasons in the face of global change: phenological research in the 21st century. *Front. Ecol. Environ.* **7**, 253–260 (2009).

26. Miralles, D. G., Gentine, P., Seneviratne, S. I. & Teuling, A. J. Land-atmospheric feedbacks during droughts and heatwaves: state of the science and current challenges. *Ann. N. Y. Acad. Sci.* **1436**, 19–35 (2019).

27. Zhang, Y. et al. Multi-decadal trends in global terrestrial evapotranspiration and its components. *Sci. Rep.* **6**, 19124 (2016).

28. Cheng, L. et al. Recent increases in terrestrial carbon uptake at little cost to the water cycle. *Nat. Commun.* **8**, 110 (2017).

29. Zeng, Z., Peng, L. & Piao, S. Response of terrestrial evapotranspiration to Earth's greening. *Curr. Opin. Environ. Sustain.* **3**, 9–25 (2018).

30. Piao, S. et al. Plant phenology and global climate change: current progresses and challenges. *Glob. Change Biol.* **25**, 1922–1940 (2019).

31. Forzieri, G., Alkama, R., Miralles, D. G. & Cescatti, A. Response to Comment on "Satellites reveal contrasting responses of regional climate to the widespread greening of Earth". *Science* **360**, eaap9664 (2018).

32. Mokany, K., Raison, R. J. & Prokushkin, A. S. Critical analysis of root: shoot ratios in terrestrial biomes. *Glob. Change Biol.* **12**, 84–96 (2006).

33. Chen, C. et al. China and India lead in greening of the world through land-use management. *Nat. Sustain.* **2**, 122–129 (2019).

34. Le Quéré, C. et al. Global Carbon Budget 2018. *Earth Syst. Sci. Data* **10**, 2141–2194 (2018).

35. Lian, X. et al. Partitioning global land evapotranspiration using CMIP5 models constrained by observations. *Nat. Clim. Change* **8**, 640–646 (2018).

36. Buermann, W. et al. Widespread seasonal compensation effects of spring warming on northern plant productivity. *Nature* **562**, 110–114 (2018).

37. Jiménez, C. et al. Global intercomparison of 12 land surface heat flux estimates. *J. Geophys. Res. Atmos.* **116**, D02102 (2011).

38. Mahowald, N. et al. Projections of leaf area index in earth system models. *Earth Syst. Dynam.* **7**, 211–229 (2016).

39. Dai, A. Increasing drought under global warming in observations and models. *Nat. Clim. Change* **3**, 52–58 (2013).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© Crown/European Union 2020
## Methods

**Vegetation dynamics.** Three satellite-based LAI products were used to analyse the changes in global vegetation for the period 1982–2016 derived from the Global Land Cover Facility (GLASS v.3, http://ftp.glcf.umd.edu/)¹⁰, the Global Inventory Modeling and Mapping Studies Normalized Difference Vegetation Index (GIMMS3g v.1, http://sites.bu.edu/cliveg/datacodes/)¹¹ and the NOAA Climate Data Record (TCDR v.4, https://www.ncei.noaa.gov/data/avhrr-land-leaf-area-index-and-fapar/access/)¹². The monthly mean 1° LAI was calculated by averaging LAI values from the original spatiotemporal resolution of each product. The residual data gaps were filled using the Harmonic Analysis of Time Series (HANTS) method¹³. The growing-season-averaged LAI was used as a proxy of vegetation growth in this study. To this end, a climatological growing season that spanned months with at least 75% of days in the greenness phase was derived from the Vegetation Index and Phenology (VIP) satellite-based product¹³ (https://vip.arizona.edu/vipdata/V4/DATAPOOL/PHENOLOGY/) and used as a reference period to derive a multi-year time series of growing-season LAI. Nonparametric Mann–Kendall trend tests were then computed at the pixel level after averaging the LAI value across a 3° × 3° spatial moving window separately for each product (Supplementary Fig. 3). The moving window aimed to preserve the spatial consistency with the retrievals of sensitivity of energy partitioning terms to LAI (see next section). A sensitivity analysis of interannual variation and the trend in LAI for the different thresholds used to identify the growing season was performed to corroborate the robustness of our results (Supplementary Fig. 4).

Results are explored for different climate zones, derived from the Köppen–Geiger world map of climate classification¹⁴, and for vegetated types (V) including forests (broad-leaved and needle-leaved trees), natural grasses and croplands. We used datasets of PFTs derived from the annual land-cover maps of the European Space Agency's Climate Change Initiative (ESA-CCI, https://www.esa-landcover-cci.org/)¹⁵ over the 2000–2014 period, referring to a simplified aggregation scheme based on physiognomy alone. On the basis of these classifications, we derived the multi-annual average cover fraction of each vegetated class (F⁶). Desert and semi-desert areas with an average growing-season LAI <0.15 m² m⁻² were excluded from the analyses.

**Energy partitioning terms and evapotranspiration products.** We focused the analysis on the interplay between interannual variations in LAI and the terms of the surface energy partitioning, including LE, H and B over the 1982–2016 period. LE was derived by combining latent heat of vaporization and ET estimates derived from four different observation-based datasets including the Global Land Evaporation Amsterdam Model (GLEAM v.3.2a, https://www.gleam.eu/)¹⁶,¹⁷, the Model Tree Ensemble (MTE, https://www.bgc-jena.mpg.de/geodb/projects/Home.php)¹⁸, the Process-based Land Surface Evapotranspiration/Heat Fluxes (PLSH, http://files.ntsg.umtd.edu/data/)¹⁹ and the Breathing Earth System Simulator (BESS, http://environment.snu.ac.kr/bess_flux/)²⁰. The BESS model represents a new long-term ET product that was specifically developed for this study and retrieved from a consolidated process-based model in combination with the three LAI satellite products that were used here. H was obtained from the closure of the energy balance by subtracting LE from the surface net radiation (RN). RN was retrieved from the ERA-interim reanalysis data²¹ (http://apps.ecmwf.int/datasets/). As such, the estimate of sensible heat implicitly includes the heat storage in canopy air and biomass and ground heat terms of the energy-balance equation. Transpiration and evapotranspiration data generated from GLEAM were also used to explore the trend in the fraction of transpired water to the total evapotranspiration (Tr/ET) (Extended Data Fig. 3i,j). Sensible fluxes derived from MTE²⁵ were used to verify the consistency of H estimates derived from the closure of the energy balance.

CO₂ concentrations may have a role in the interplay between changes in LAI and energy partitioning by leading to the partial closure of stomata and restricting the diffusion of water vapour out of leaves¹⁵. The use of PLSH and BESS, which explicitly integrate the CO₂ effect on stomatal conductance, enables us to account for such mechanisms. On the other hand, LE estimates based on GLEAM are independent from LAI (or other optical remote-sensing metrics), which preserves the interplay between LE and LAI estimates from possible circularity effects. The vegetation status in GLEAM is characterized by the vegetation optical depth, a microwave-based vegetation parameter related to vegetation water content and biomass. By contrast, MTE, PLSH and BESS are based on a set of satellite-based predictors, including the normalized difference vegetation index, the fraction of photosynthetically active radiation and LAI (Supplementary Table 2).

**Inferring supply and demand limitations of LE.** We inferred the primary limitations of LE, atmospheric demand or moisture supply, by comparing the correlations between the growing-season averages of LE and growing-season averages of temperature and precipitation, consistent with previous studies¹⁴,³². Growing-season dates were derived from surface phenological data from VIP¹⁴ (as computed for LAI). As temperature, radiation and the vapour-pressure deficit are strongly correlated, temperature can be used as a proxy for atmospheric demand. We then compared the correlation between LE and precipitation (ρ(LE, P)) and the correlation between LE and temperature (ρ(LE, T))—both of which were computed

over the whole time series spanning the 1982–2016 period—and defined supply-limited zones for which ρ(LE, P) > ρ(LE, T) and demand-limited zones for which ρ(LE, P) < ρ(LE, T). We then obtained a 'static' classification map (moisture supply-limited compared with atmospheric demand-limited), resulting from the average of the multiple correlation maps obtained from the different LE (ET) products and climate data used in this study. Such clustering was used in combination with the vegetation map for which only pixels with ≥80% of the vegetated cover and ≤10% of the irrigated area were included in the study domain (Supplementary Fig. 5). We therefore used the Global Map of Irrigation Areas (GMIA, http://www.fao.org/irn/water/aquastat/irrigationmap/index.stm) derived from statistical census data⁵³ for the year 2005. Furthermore, to explore the temporal evolution of supply- and demand-limited zones, we replicated the above-mentioned classification scheme across moving windows of different length (7 and 13 years). These analyses produced a set of 'dynamic' classification maps (Extended Data Fig. 4).

**Sensitivity of surface energy partitioning on LAI changes.** The sensitivity of surface energy partitioning terms (LE, H or B; hereafter referred to as Z for simplicity) to changes in LAI was initially computed at the annual scale for the whole 1982–2016 period. It was expressed as the partial derivative that results from a multiple linear regression, relating the interannual differences in the Z component (ΔZ) to interannual differences in the growing-season-averaged LAI (ΔLAI), annually averaged air temperature (ΔT), annually cumulated precipitation (ΔP) and annually averaged incoming shortwave radiation (ΔSWin):

$$ \Delta Z = \beta_0 + \frac{\partial Z}{\partial LAI} \Delta LAI + \frac{\partial Z}{\partial T} \Delta T + \frac{\partial Z}{\partial P} \Delta P + \frac{\partial Z}{\partial SW_{in}} \Delta SW_{in} \quad (1) $$

T, SWin and P are the climate data used to derive each ET product (Supplementary Table 2), whereas LAI is retrieved from the satellite products described in the previous section (GLASS v.3, GIMMS3g v.3 and TCDR v.4). All data have been linearly resampled before analysis at the common 1° × 1° spatial resolution. We applied equation (1) to each unique combination of energy fluxes and the LAI dataset, therefore resulting in a 12-member ensemble of sensitivity estimates for each Z term. Such an approach based on the difference between two consecutive years (the Δ operator) disentangles the resulting signal from possible long-term dependencies on covariates (for example, the combined effect of rising temperatures and CO₂ concentrations on long-term LAI trends, as well as effects of long-term drying on soil moisture and biomass). The derived signal $\frac{\partial Z}{\partial LAI}$ integrates the bidirectional interactions between LAI and the Z term. To better sample the statistical inferences all predictors in equation (1) were quantified for each pixel over a centred 3° × 3° spatial window. This approach factors out the effects of LAI on energy fluxes that are triggered by a variation in the main climate drivers (for example, increase in P triggering larger LAI and LE fluxes). However, the regression model in equation (1) assumes a linear interplay between the response variable and predictors and does not account for the possible covariation among predictors. The use of this method, instead of more sophisticated techniques, appears a reasonable approach considering: (1) the length of the time series (35 years); (2) the choice of a parsimonious approach that can be applied consistently across different products and variables; (3) the ability to capture possible emerging first-order temporal changes in the signal.

We also applied equation (1) to monthly timescales separately for each LAI product. The sensitivity at the monthly timescale was computed by using growing-season averaged LAI values and monthly scale climate drivers to minimize the potential biases of satellite retrieval of LAI in snow-cover conditions and to explore the effects of the changes in background climate. Note that monthly scale Δ values for climate drivers are calculated as the difference between the same months of two consecutive years. Furthermore, to better characterize soil moisture conditions, monthly scale P accounts for concurrent and lagged cumulated precipitation for which the contributions are derived from an empirically derived decay exponential function assuming a soil depth¹⁴ of 1 m as follows:

$$ P_{t_0} = \sum_{i=0}^{n} \frac{P_i}{t_i - t_{i+1}} \left(1 - e^{-(t_i - t_{i+1})}\right) e^{-(t_i - t_0)} \quad (2) $$

where t₀ and tᵢ refer to the current month and the ith lagged month, respectively.

**Temporal variations in the sensitivity of energy partitioning terms to LAI changes.** To explore possible long-term variations in the sensitivity of energy partitioning terms to LAI changes, equation (1) was also computed on annual scale over 7-year and 13-year temporal moving windows. This analysis was complemented with sensitivities estimated for two consecutive independent periods ranging from 1982 to 1999 (tᵢ) and from 2000 to 2016 (t₂). Absolute (Δabs) and relative (Δrel) changes in sensitivities were quantified and a t-test was then used to determine whether the two samples $\frac{\partial Z}{\partial LAI}$ and $\frac{\partial Z}{\partial LAI}$, were significantly different from each other. Temporal variations in the sensitivity of sensible fluxes to LAI were also computed by using native estimates (Hₙ) directly provided by MTE. Results of this latter comparison are largely consistent with estimates obtained from H = RN – LE, confirming the marginal effects of residual heat-storage flux (Supplementary Fig. 6).
**Potential effects of changes in satellite sensors.** The LAI datasets used in this study have been generated from time series of satellite observations that have been specifically harmonized to remove any biases caused by changes in the sensors. However, we cannot exclude that some residual effects of sensor change could still influence the year-to-year variations in LAI and thus the sensitivity of energy partitioning to LAI changes<sup>23,35</sup>. We therefore explored the variations in sensitivity recorded during the two aforementioned periods, 1982–1999 (t<sub>1</sub>) and 2000–2016 (t<sub>2</sub>) separately for each LAI product (Supplementary Fig. 1). The split between these two periods reflects a major change in the monitoring system given the first year of operation of Terra MODIS in 2000—a sensor used in several LAI products in the t<sub>2</sub> period. During the t<sub>1</sub> period, LAI data were exclusively based on AVHRR acquisitions.

We evaluated the potential effects of changes in sensors by exploring the climate control on sensitivity ($\partial Z/\partial LAI$) for the t<sub>1</sub> and t<sub>2</sub> periods. Therefore, we expressed the sensitivity estimates as a function of the aridity index (AI)<sup>56</sup>, which was quantified as:

$$AI = \frac{P}{T + 33} \qquad (3)$$

where P and T are the climatological estimates of annual cumulated precipitation and annual average temperature computed for the reference temporal period (Extended Data Fig. 2a,c,e). Based on the Kolmogorov–Smirnov significance test (P-value < 0.05), we cannot reject the hypothesis that the two resulting samples come from the same distribution. In case of a systematic bias on the climatic control on $\partial Z/\partial LAI$ the two curves should have been statistically different.

Furthermore, we derived the sensitivity of energy fluxes to LAI changes for the whole 1982–2016 period using a spatial gradient derived from a climatic space instead of temporal variability. The sensitivities were therefore binned in a precipitation–temperature (PT) space in which every bin was equally affected by variations in the sensors, and for which the gradients were therefore independent of sensor changes over time (Fig. 1a,d,g). We then extrapolated annual sensitivity values ($\partial Z/\partial LAI$) from the PT domain based on the annual precipitation and temperature for each grid of the globe and for each year. Therefore, the temporal changes in sensitivity derived with this second methodology are fully conditioned on the changes in climate (Extended Data Fig. 2b,d,f). The emerging temporal changes in $\partial Z/\partial LAI$ show an increase in sensitivity of energy partitioning to LAI consistent with $\partial Z/\partial LAI$, thus corroborating the role of environmental conditions as drivers of the increase in vegetation control on energy partitioning.

**Effects of LAI trends on available surface energy partitioning.** Variations in surface energy partitioning associated with the long-term variations in LAI ($\delta Z^{LAI}$) have been computed by applying the previously described methodology<sup>12</sup> and are expressed by the following equation:

$$\delta Z^{LAI} = \frac{\partial Z}{\partial LAI} \delta LAI \qquad (4)$$

where $\delta LAI$ is the long-term trend in growing-season-averaged LAI and $\partial Z/\partial LAI$ is the sensitivity of Z to LAI (equation (1)) quantified over the 1982–2016 period (both at annual and monthly levels). Consistent with the assessment of sensitivity, a 12-member ensemble of trend estimates was derived based on the different combinations of original ET and LAI products of each energy term. We implicitly assumed that the sensitivity $\partial Z/\partial LAI$ computed at the interannual scale (equation (1)) is an appropriate metric to estimate the net climate impact on the phenomena. However, at longer timescales (from decades to a century) additional ecosystem processes may emerge, such as adaptation phenomena driven by the change in species and shifting biomes, which could affect the future trends of the sensitivity. An indication of the long-term sensitivity in case of full adaptation can be derived from the analysis of the sensitivity in the spatial domain (Fig. 1a,d,g and Extended Data Fig. 2b,d,f).

To compare the effects of LAI on the surface energy partitioning terms with those resulting from other drivers of the regression, equation (4) was similarly applied to the sensitivities and long-term trend estimates of annually averaged T, P and SW<sub>in</sub>. The marginal effects of each single predictor have been quantified and compared with the overall trends in surface energy fluxes. The fraction of the overall trend in the energy term explained by a given predictor is then quantified as the ratio between the predictor-specific trend and the overall trend, which were both averaged globally.

**Disentangling human land-use management.** We disentangled the marginal contribution of forests, natural grasses and croplands to the global signal of the long-term LAI effect on energy partitioning ($\delta Z^{LAI}$) computed at the annual scale for the whole 1982–2016 period. The marginal contribution of a given vegetation type ($V^V$) is derived by using the grid-cell cover fractions ($F^V$) as weights, based on the following equation:

$$Q^V (\delta Z^{LAI}) = 100 \frac{\sum_{i=1}^{n} (F_i^V A_i \delta Z^{LAI}_i)}{\sum_{i=1}^{n} A_i \delta Z^{LAI}_i} \qquad (5)$$

where i represents a pixel, n is the total number of vegetated pixels globally, $\delta Z^{LAI}$ is the effect of a pixel, A<sub>i</sub> is the area of a pixel that varies with latitudes. A similar

approach was used to derive the marginal contribution of the different vegetation types on the global land area, the long-term trend in growing-season-averaged LAI (δLAI) and the sensitivity of the energy partitioning terms to LAI changes ($\partial Z/\partial LAI$). To characterize the biome-specific modulation effects on the LAI–energy interplay, we binned the cover fractions of each vegetation type as a function of δLAI and $\partial Z/\partial LAI$.

Furthermore, we derived a set of different land-use management conditions by excluding all pixels with a fraction of croplands that varied from 30% to 70% from the study. We analysed the relative changes in sensitivities over time ($\Delta rel$) and the climate control on the effects ($\delta Z^{LAI}$) for different land-use management conditions (cropland percentage). The climate control was determined by binning $\delta Z^{LAI}$ as a function of the aridity index (equation (3)) and by testing the similarity of curves by the Kolmogorov–Smirnov test. Spatial and temporal patterns of sensitivity ($\partial Z/\partial LAI$) and effects ($\delta Z^{LAI}$) appear to be substantially independent of the crop coverage (Supplementary Fig. 7), thus confirming the climate controls on the LAI–energy interplay even in vegetated lands that are subject to human land-use management.

**Simulations of the LSMs.** To complement the analysis based on observational products, we used simulations from ten state-of-the-art LSMs completed within the TRENDY v7 project<sup>34</sup>, including CABLE-POP, CLASS-CTEM, CLM5.0, DLEM, ISAM, JSBACH, JULES, LPX, ORCHIDEE-CNP and VISIT. All models provide prognostic estimates of LAI and LE (only two models provided H). To analyse the modelled relative contributions of external factors to changes in the sensitivity of LE to LAI, we used factorial simulations obtained for the 1982–2016 period under different scenarios: changes in CO<sub>2</sub>, climate and land use (S3, the most realistic scenario); changes in CO<sub>2</sub> only (S1) and changes in climate and land use only (S3–S1). For each run, we quantified the modelled sensitivity of LE to changes in LAI (equation (1)) and the associated long-term effect (equation (4)). Results for the climate and land-use change scenario are obtained by subtracting the sensitivity computed under S1 to sensitivity computed under S3. Model results are compared with analogous estimates derived from satellite-observation-based products and bias patterns explored across the gradient of the differences in absolute ΔLAI between data and models.

**Multi-product ensembles.** To better capture the emerging signals and account for the possible differences across ET and LAI products, we calculated the multi-product average trend in the surface energy partitioning Z term of the 12 experiments obtained for different observation-driven LAI–ET combinations. To derive a global estimate of the trend and related uncertainty and to fulfil the assumption of uncorrelated errors, we subset the global domain-extracting pixels with non-overlapping spatial windows—that is, only 3 equally spaced pixels for the latitude and 3 pixels for the longitude were selected. The sampling was replicated nine times to progressively cover the full global domain. For each of the nine global subsets (D), we derived the zonal median of the trend, weighting each grid-cell value based on its area. The global estimates of the trend in the Z term were quantified as the average of the estimates derived from the ensemble of the nine global subsets. The analysis was replicated separately for each experiment. We then derived the average and the corresponding s.e. of the ensemble of single-experiment global estimates. Similar procedures were used to quantify the multi-product ensemble average and uncertainty of sensitivities and trend in the growing-season-averaged LAI. The same approach was used for the ensemble of LSMs. We refer to the average and s.e. of the ensembles in the text and figures where not differently indicated.

**Reporting Summary.** Further information on research design is available in the Nature Research Reporting Summary linked to this article.

## Data availability

The observation-driven datasets analysed in this study are publicly available as referenced within the Article. Simulations from ten LSMs (CABLE-POP, CLASS-CTEM, CLM5.0, DLEM, ISAM, JSBACH, JULES, LPX-Bern, ORCHIDEE-CNP and VISIT) are available from the TRENDY dataset on request from S.S. All generated data are available from the corresponding author on request.

## Code availability

The custom MATLAB (R2017b) code written to read and analyse data and generate figures is fully available on request from the corresponding author.

## References

40. Xiao, Z. et al. Long-time-series global land surface satellite leaf area index product derived from MODIS and AVHRR surface reflectance. *IEEE Trans. Geosci. Remote Sens.* **54**, 5301–5318 (2016).

41. Zhu, Z. et al. Global data sets of vegetation leaf area index (LAI)3g and fraction of photosynthetically active radiation (FPAR)3g derived from global inventory modeling and mapping studies (GIMMS) normalized difference vegetation index (NDVI3g) for the period 1981 to 2011. *Remote Sens.* **5**, 927–948 (2013).
42. Claverie, M., Matthews, J. L., Vermote, E. F. & Justice, C. O. A 30+ Year AVHRR LAI and FAPAR Climate data record: algorithm description and validation. *Remote Sens.* **8**, 263 (2016).

43. Zhang, X., Tan, B. & Yu, Y. Interannual variations and trends in global land surface phenology derived from enhanced vegetation index during 1982–2010. *Int. J. Biometeorol.* **58**, 547–564 (2014).

44. Kottek, M., Grieser, J., Beck, C., Rudolf, B. & Rubel, F. World map of the Köppen–Geiger climate classification updated. *Meteorol. Z.* **15**, 259–263 (2006).

45. ESA. *Land CoverCCI Product User Guide* v.2 (2017).

46. Miralles, D. G. et al. Global land-surface evaporation estimated from satellite-based observations. *Hydrol. Earth Syst. Sci.* **15**, 453–469 (2011).

47. Martens, B. et al. GLEAM v3: satellite-based land evaporation and root-zone soil moisture. *Geosci. Model Dev.* **10**, 1903–1925 (2017).

48. Jung, M. et al. Global patterns of land-atmosphere fluxes of carbon dioxide, latent heat, and sensible heat derived from eddy covariance, satellite, and meteorological observations. *J. Geophys. Res. Biogeosci.* **116**, G00J07 (2011).

49. Zhang, K. et al. Vegetation greening and climate change promote multi-decadal rises of global land evapotranspiration. *Sci. Rep.* **5**, 15956 (2015).

50. Jiang, C. & Ryu, Y. Multi-scale evaluation of global gross primary productivity and evapotranspiration products derived from Breathing Earth System Simulator (BESS). *Remote Sens. Environ.* **186**, 528–547 (2016).

51. Dee, D. P. et al. The ERA-Interim reanalysis: configuration and performance of the data assimilation system. *Q. J. R. Meteorol. Soc.* **137**, 553–597 (2011).

52. Jung, M. et al. Recent decline in the global land evapotranspiration trend due to limited moisture supply. *Nature* **467**, 951–954 (2010).

53. Siebert, S. et al. A global data set of the extent of irrigated land from 1900 to 2005. *Hydrol. Earth Syst. Sci.* **19**, 1521–1545 (2015).

54. Pan, F., Peters-Lidard, C. D. & Sale, M. J. An analytical method for predicting surface soil moisture from rainfall observations. *Water Resour. Res.* **39**, 1314 (2003).

55. Tian, F. et al. Evaluating temporal consistency of long-term global NDVI datasets for trend analysis. *Remote Sens. Environ.* **163**, 326–340 (2015).

56. Köppen, W. *Die Klimate der Erde: Grundriss der Klimakunde* (Walter de Gruyter und Co, 1923).

## Acknowledgements

The study was funded by the FP7 LUC4C project (grant number 603542). D.G.M. acknowledges funding from the European Research Council (ERC) under grant agreement 715254 (DRY-2-DRY). P.C. acknowledges support from European Research Council Synergy project SyG-2013-610028 IMBALANCE-P and ANR (reference ANR-16-CONV-0003 (CLAND)).

## Author contributions

G.F. and A.C. conceived and designed the study; D.G.M. and B.M. provided GLEAM data; C.J. and Y.R. produced the archive of long-term BESS data and harmonized LAI datasets; K.Z. provided PLSH data; A.W., A.A., D.S.G., V.K.A., S.L., D.L., E.K., J.E.M.S.N., H.T., P.F. and S.S. ran the TRENDY v.7 simulations; R.A. harmonized LSM simulations; G.F. analysed the data, G.F. and A.C. interpreted the results and wrote the manuscript with contributions from all co-authors.

## Competing interests

The authors declare no competing interests.

## Additional information

*Extended data* is available for this paper at https://doi.org/10.1038/s41558-020-0717-0.

*Supplementary information* is available for this paper at https://doi.org/10.1038/s41558-020-0717-0.

*Correspondence and requests for materials* should be addressed to G.F.

*Peer review information* *Nature Climate Change* thanks Liang Chen and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

*Reprints and permissions information* is available at www.nature.com/reprints.
<figure><img src="image_32.png" /></figure>

**Extended Data Fig. 1** | Temporal variations of sensitivity of LE to LAI changes for single LAI and ET products and ensemble averages. Temporal variations of sensitivities computed over a 13-year moving window for supply- and demand limited regions and the whole globe. Labels report the relative changes in sensitivities (Δrel, Methods) between the 1982-1999 period and the 2000-2016 period, '*' indicates the Student's t-test significance with p-value ≤ 0.05. (a-d, f-i, k-n) represent single combinations of LAI and ET products. (p-s) represent ensemble results averaged across multiple LAI products. (e, j, o) represent ensemble results averaged across multiple ET products. (t) represents ensemble results across multiple LAI and ET products. In the ensemble configurations, values are expressed as ensemble average ± standard error. The spatial domains of supply- and demand-limited regions are shown in Supplementary Fig. 5.
<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>

**Extended Data Fig. 2** | Climate control on the vegetation-energy interplay. (a, c, e) Sensitivity of latent heat ($\frac{\partial LE}{\partial LAI}$), sensible heat ($\frac{\partial H}{\partial LAI}$) and Bowen ratio ($\frac{\partial B}{\partial LAI}$) to LAI changes (on the y-axis) separately computed for the 1982-1999 and 2000-2016 periods and binned as a function of the aridity index (on the x-axis) (Eq. (3), Methods). Results of the Kolmogorov-Smirnov test is shown in label and reflects the significance level ($p_{ks}$) to reject the null hypothesis of dissimilar curves. (b, d, f) Temporal variations of sensitivities extrapolated from changes in precipitation and climate for moisture supply- and atmospheric demand-limited regions and the whole globe displayed with respect to the first year sensitivity (year 1982). Numbers at the bottom of the panel report the sensitivity values of the first year, while numbers on top report the coefficient of regression (R²) and the significance (Mann-Kendall test; $p_{mk}$) of the fitting linear regression models. The spatial domains of supply- and demand-limited regions are shown in Supplementary Fig. 5.
<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /></figure>

<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>

<figure><img src="image_44.png" /></figure>

<figure><img src="image_45.png" /></figure>

<figure><img src="image_46.png" /></figure>

<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /></figure>
**Extended Data Fig. 3** | Trends in vegetation and climate drivers. Spatial patterns (**a**) and climate space (**b**) of long-term trend (1982–2016) in growing season averaged leaf area index ($\delta$LAI). (**c**, **d**), (**e**, **f**), (**g**, **h**) and (**i**, **j**) as (**a**, **b**) but for temperature ($\delta$T), precipitation ($\delta$P), short-wave incoming radiation ($\delta$SW<sub>IN</sub>) and the ratio between transpiration and total ET ($\delta$Tr/ET). Values reflect the ensemble average of multiple products used in this study (e.g., GLASS v3, GIMMS3g v3 and TCDR v4 for $\delta$LAI). Areas in (**a**, **c**, **e**, **g**, **i**) labelled with black dots indicate trends that are statistically significant (Mann-Kendall test; p-value$\le$0.05). Values in (**b**, **d**, **f**, **h**, **j**) are binned as a function of climatological mean precipitation (*P*, on the x-axis) and air temperature (*T*, on the y-axis) and black dots show bins with average values statistically different from zero (Student's t-test; p-value$\le$0.05).
<figure><img src="image_49.png" /></figure>

<figure><img src="image_50.png" /></figure>

**Extended Data Fig. 4** | Dynamics of global coverage of moisture supply- and atmospheric demand-limited regions. (a) Temporal variations in global fraction of moisture supply- and atmospheric demand-limited zones, shown respectively on the left and right y-axis, based on the ensemble average of observation-driven products (Methods). (b) as (a) but for the ensemble average of LSM simulations.
<figure><img src="image_51.png" /></figure>

**Extended Data Fig. 5** | Sensitivity of energy partitioning terms to multiple drivers. (a, d, g, j) Sensitivity of latent heat to changes in LAI ($\frac{\partial LE}{\partial LAI}$), temperature ($\frac{\partial LE}{\partial T}$), precipitation ($\frac{\partial LE}{\partial P}$) and short-wave incoming radiation ($\frac{\partial LE}{\partial SW_{IN}}$) calculated for the 1982–2016 period as ensemble average of all observation-driven estimates. (b, e, h, k) and (c, f, i, l) as (a, d, g, j) but for sensible heat (H) and Bowen ratio (B). Areas labelled with black dots indicate estimates that are statistically significant (Student's t-test; p-value ≤0.05).
<figure><img src="image_52.png" /><figcaption>Extended Data Fig. 6 | See next page for caption.</figcaption></figure>
**Extended Data Fig. 6** | Temporal variations of sensitivity of LE to LAI changes for single land surface models and ensemble averages. Temporal variations of sensitivities computed over a 13-year moving window for supply- and demand-limited regions and the whole globe displayed with respect to the first year sensitivity (year 1989). Results refer to different scenarios: changes in CO₂, climate and land use (S3, the most realistic scenario, left column); changes in CO₂ only (S1, middle column) and changes in climate and land use only (S3-S1, right column). Black labels report the relative changes (Δrel, Methods) in global sensitivities between the 1982-1999 period and the 2000-2016 period, '*' indicates the Student's t-test significance with p-value≤0.05, while the number in brackets refer to the first year sensitivity value. Results are shown for each single land surface model (a1-a10, b1-b10, c1-c10) and for the ensemble average (a11, b11, c11). The spatial domains of supply- and demand-limited regions are shown in Supplementary Fig. 5.
<figure><img src="image_53.png" /><figcaption>Extended Data Fig. 7 | See next page for caption.</figcaption></figure>
**Extended Data Fig. 7** | Latitudinal profiles of sensitivity of LE to LAI, long-term trend in LAI, and the resulting greening effect on LE for each land surface model and the ensemble average. Latitudinal profiles of sensitivity of latent heat to changes in LAI (δLAI / δLAI) (left column), long-term trends in LAI (δLAI) (middle column) and associated changes in long-term trend in LE (δLE/LAI) (right column). Results are shown for each single land surface model (a1-a10, b1-b10, c1-c10) and for the ensemble average (a11, b11, c11). In the last row, results from observation-driven products are shown for comparison (red lines). Background colours in (c11) show areas of the globe subject to compensatory effects (grey), additive positive biases (orange) and additive negative biases (blue).
# Reporting Summary

Nature Research wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Research policies, see *Authors & Referees* and the *Editorial Policy Checklist*.

## Statistics

For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a | Confirmed

☐ ☑ The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement

☐ ☑ A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly

☐ ☑ The statistical test(s) used AND whether they are one- or two-sided
*Only common tests should be described solely by name; describe more complex techniques in the Methods section.*

☐ ☑ A description of all covariates tested

☐ ☑ A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons

☐ ☑ A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals)

☐ ☑ For null hypothesis testing, the test statistic (e.g. *F*, *t*, *r*) with confidence intervals, effect sizes, degrees of freedom and *P* value noted
*Give *P* values as exact values whenever suitable.*

☒ ☐ For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings

☒ ☐ For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes

☐ ☑ Estimates of effect sizes (e.g. Cohen's *d*, Pearson's *r*), indicating how they were calculated

*Our web collection on statistics for biologists contains articles on many of the points above.*

## Software and code

Policy information about [availability of computer code](#)

Data collection
Phyton was used to harmonize simulations of land surface models
Matlab was used to harmonize satellite data

Data analysis
Matlab was used to perform all the analyses described in the manuscript

For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors/reviewers.
We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research guidelines for submitting code & software for further information.

## Data

Policy information about [availability of data](#)

All manuscripts must include a [data availability statement](#). This statement should provide the following information, where applicable:

- Accession codes, unique identifiers, or web links for publicly available datasets

- A list of figures that have associated raw data

- A description of any restrictions on data availability

The observation-driven datasets analyzed in this study are publicly available as referenced within the article. Simulations from ten Land Surface Models (CABLE-POP, CLASS-CTEM, CLM5.0, DLEM, ISAM, JSBACH, JULES, LPX-Bern, ORCHIDEE-CNP, VISIT) are available from the TRENDY dataset via a request to S. Sitch. All generated data and scripts are available from the corresponding author on request.
## Field-specific reporting

Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection.

<figure><img src="image_54.png" /></figure>

Life sciences

<figure><img src="image_55.png" /></figure>

Behavioural & social sciences

<figure><img src="image_56.png" /></figure>

Ecological, evolutionary & environmental sciences

For a reference copy of the document with all sections, see [nature.com/documents/nr-reporting-summary-flat.pdf](https://nature.com/documents/nr-reporting-summary-flat.pdf)

## Ecological, evolutionary & environmental sciences study design

All studies must disclose on these points even when the disclosure is negative.

Study description

Our analysis provides observational evidence that changes in vegetation density during the past three and half decades (1982-2016) have affected substantially the terrestrial energy fluxes, by enhancing latent heat over sensible heat, ultimately leading to an intensification of the terrestrial water cycle. Notably, the sensitivity of evapotranspiration to variation in leaf area index has sharply increased in the observation period, suggesting an enhanced biological control over the terrestrial water cycle.

Research sample

Our analysis focuses on supply- and demand limited regions and the whole globe in order to elucidate the key physical mechanisms driving the temporal changes in vegetation control on surface energy partitioning. We also distinguish the vegetation response between grass and tree biomes to investigate the potential dependence on vegetation structure.

Sampling strategy

Our analysis explores the vegetation control on surface energy partitioning at 1 degree spatial resolution and at annual and monthly scales. These spatial and temporal scales are appropriate to quantify the short-term (interannual) changes in biophysical processes considered in this study.

Data collection

G.F. (first and corresponding author) collected satellite data and land surface model output. Satellite data were retrieved from open-access archives and land surface model simulations directly from the modeling teams, as reported in the data availability paragraph.

Timing and spatial scale

Our analysis explores the vegetation control on surface energy partitioning at global scale over the 1982-2016 period.

Data exclusions

We focus on vegetated pixels with cover fraction ≥80%, growing season averaged LAI≥ 0.15 m2m-2 and with ≤10% of irrigated. Pixels out of this domain are excluded because the signal in these regions is not robust and may bias the final results.

Reproducibility

All methodological procedures employed in the study are described extensively in the Method section to allow reproducibility of our results (multiple linear regression models, spatial aggregation, error propagation, interannual variations).

Randomization

Randomization is not relevant for this study because we provide estimates based on an ensemble of satellite products and land surface models developed and forced with reasonably independent approaches and datasets.

Blinding

Blinding was not relevant for this study because we explore the whole globe based on satellite observations and land surface models

Did the study involve field work?

<figure><img src="image_57.png" /></figure>

Yes

<figure><img src="image_58.png" /></figure>

No

## Reporting for specific materials, systems and methods

We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response.

### Materials & experimental systems

<table><thead><tr><th>n/a</th><th>Involved in the study</th></tr></thead><tbody><tr><td>☒</td><td>Antibodies</td></tr><tr><td>☒</td><td>Eukaryotic cell lines</td></tr><tr><td>☒</td><td>Palaeontology</td></tr><tr><td>☒</td><td>Animals and other organisms</td></tr><tr><td>☒</td><td>Human research participants</td></tr><tr><td>☒</td><td>Clinical data</td></tr></tbody></table>

### Methods

<table><thead><tr><th>n/a</th><th>Involved in the study</th></tr></thead><tbody><tr><td>☒</td><td>ChIP-seq</td></tr><tr><td>☒</td><td>Flow cytometry</td></tr><tr><td>☒</td><td>MRI-based neuroimaging</td></tr></tbody></table>