## Abstract

Satellite observations show widespread increasing trends of leaf area index (LAI), known as the Earth greening. However, the biophysical impacts of this greening on land surface temperature (LST) remain unclear. Here, we quantify the biophysical impacts of Earth greening on LST from 2000 to 2014 and disentangle the contributions of different factors using a physically based attribution model. We find that 93% of the global vegetated area shows negative sensitivity of LST to LAI increase at the annual scale, especially for semiarid woody vegetation. Further considering the LAI trends (*P* ≤ 0.1), 30% of the global vegetated area is cooled by these trends and 5% is warmed. Aerodynamic resistance is the dominant factor in controlling Earth greening’s biophysical impacts: The increase in LAI produces a decrease in aerodynamic resistance, thereby favoring increased turbulent heat transfer between the land and the atmosphere, especially latent heat flux.

## INTRODUCTION

Vegetation is a key regulator of land-atmosphere exchanges of heat, mass, and momentum (*1*, *2*). Satellite records of leaf area index (LAI), a measure of the aboveground vegetation abundance, indicate widespread increasing trends since the 1980s (*3*), due to warming (*4*, *5*), CO_{2} fertilization (*6*), and land management (*7*, *8*). The Earth greening directly affects the surface energy budget through influencing various biophysical factors: (i) albedo (α), which controls the fraction of solar radiation absorbed by the surface; (ii) aerodynamic resistance (*r*_{a}), which characterizes the efficiency of turbulent transfer of heat and water vapor; (iii) surface resistance (*r*_{s}), which is the additional resistance to water vapor transport through the soil and the pores on leaves; and (iv) emissivity (ε), which is the effectiveness of a surface in emitting and absorbing longwave radiation. Changes in these biophysical factors can strongly affect the radiometric land surface temperature (LST).

Despite their importance, the magnitude of these biophysical changes associated with Earth greening and their impacts on LST remain poorly understood (*2*, *9*–*12*). First, it is unclear which biophysical factor dominates the effects of Earth greening on LST. A global modeling study suggested that changes in aerodynamic resistance caused by the increasing LAI play a negligible role (*9*). On the other hand, several other studies demonstrate that aerodynamic resistance is the most critical biophysical factor that controls the LST response to vegetation changes (*13*–*15*). Quantifying the biophysical impacts of Earth greening on LST and addressing the dominant biophysical factor frame the scope of our study.

Second, the biophysical impacts of Earth greening on LST cannot be simply quantified by the observed changes in LST from satellites. The observed changes in LST are compound effects of Earth greening and large-scale climate change (fig. S1)—here, the Earth greening means an increasing trend in LAI, and the large-scale climate change refers to changing atmospheric conditions such as rising air temperature. Consequently, the sensitivity of LST to LAI inferred from statistical methods, e.g., multiple linear regression (*10*), may suffer causality issues [see, e.g., a recent debate on whether greening is the cause or effect of warming in the boreal zone (*10*–*12*)]. In this study, we do not aim to quantify the interplay between Earth greening and rising air temperature. Instead, we aim to quantify the impacts of Earth greening on LST independent of those from large-scale climate change.

To do so, we use a physically based attribution method, the two-resistance mechanism (TRM) method (*14*, *16*, *17*). It quantifies the sensitivities of LST to Earth greening through biophysical pathways [α, *r*_{a}, *r*_{s}, ε, and ground heat flux (*G*)] and thus the contributions of these biophysical changes to changes in LST. Unlike previous studies that compare the effects of changes in sensible and latent heat fluxes (*9*, *10*, *18*, *19*), we take one step further by considering changes in aerodynamic resistance (mostly related to surface roughness) and surface resistance (mostly related to soil moisture and vegetation characteristics), which circumvents the strong correlation between sensible and latent heat fluxes through aerodynamic resistance (*14*, *16*, *17*).

Changes in biophysical factors are caused not only by the Earth greening but also by changing atmospheric conditions. For example, changes in surface resistance could be a result of changes in LAI, air temperature, and specific humidity (*20*). Hence, to quantify the biophysical impacts of Earth greening on LST, we need to further consider the sensitivities of biophysical factors to changes in LAI (e.g., _{sat} is the satellite-derived long-term trend in LAI (fig. S2A) multiplied by the length of study period (15 years);

## RESULTS

### Sensitivity of LST to LAI through biophysical factors

The sensitivity of LST to LAI diagnosed from Community Land Model Version 5 (CLM5) outputs (hereafter ^{2} m^{−2} (mean ± 1 SD, where SD indicates spatial variability). We find that the mean magnitude of ^{2} m^{−2}) than those in high-latitude (−0.34 K m^{2} m^{−2}) and tropical regions (−0.29 K m^{2} m^{−2}) (Table 1). We further investigated ^{2} m^{−2}) and croplands (−0.43 K m^{2} m^{−2}), followed by grasslands (−0.36 K m^{2} m^{−2}), but weak in forests (−0.23 K m^{2} m^{−2}) (Table 1). Our analysis shows that OWV mostly grows in relatively dry regions where a change in LAI can substantially alter the efficiency of convective heat transfer between the land and atmosphere through changing the aerodynamic resistance (fig. S3B). A stronger sensitivity in dry regions is evidenced by the fact that the magnitude of *r*_{a}, and *r*_{s}) as well as LST where LAI is high (fig. S4, G to I) (*15*, *21*, *22*). As a result, the magnitude of ^{2} m^{−2} (−0.45 K m^{2} m^{−2}) is much larger than regions with LAI >4 m^{2} m^{−2} (−0.09 K m^{2} m^{−2}) (Table 1).

Multiple lines of evidence show that our *23*, *24*), we use the Modern-Era Retrospective analysis for Research and Applications Version 2 (MERRA-2) reanalysis data and seven fully coupled CMIP5 (the fifth phase of the Coupled Model Intercomparison Project) model outputs to compute the sensitivities of LST to biophysical factors (e.g.,

### Earth greening cools LST

The change of LST due to Earth greening (i.e., _{sat}) (Eq. 1). Our results show that *9*). On the basis of this agreement, but acknowledging that the previous study used a different model, we conjecture that the Earth greening affects LST mostly through biophysical effects with small contributions from atmospheric feedbacks. We translate ^{21} J), which is more than five times the world’s total primary energy supply in 2015 (5.71 × 10^{20} J) (*25*). From this perspective, the Earth greening–induced cooling effect is also much stronger than the warming caused by land cover change from 2000 to 2015 (1.21 × 10^{20} J) (*19*), given that one-third of the global vegetated area shows statistically significant trends in LAI (*7*).

Because the change in LST is the product of the sensitivity, _{sat} (Eq. 1), high sensitivity of LST to LAI does not necessarily imply a large change in LST, and vice versa. For example, _{sat} so that _{sat} is small (Fig. 1F and fig. S2A). In contrast, in subtropical to temperate regions of the Northern Hemisphere (eastern China, India, and Eastern Europe), _{sat} (fig. S2A).

### The dominant role of aerodynamic resistance

We quantify the contributions from different biophysical factors to *r*_{a}, *r*_{s}, ε, and *G* (Eq. 1). At the annual scale, *r*_{a} plays a dominant role in regulating the biophysical impacts of Earth greening on LST, while the impacts from ε and *G* are two to three orders of magnitude smaller than the other factors and thus negligible (*26*, *27*). In CLM5, the *r*_{a} pathway dominates 82% of the global vegetated area, followed by the *r*_{s} pathway (14.7%) and the α pathway (3.6%) (Fig. 2A). In terms of the areal fraction, 92% of the forests and 89% of OWV show that the *r*_{a} pathway plays a more prominent role than α and *r*_{s}, while the dominant fraction of *r*_{a} is reduced to 73% in grasslands and 68% in croplands (Fig. 2A). The increase in LAI results in a reduced aerodynamic resistance with a global cooling effect of −0.34 ± 0.22 K m^{2} m^{−2} (Fig. 2B), which, on average, contributes to 77% of the variance of the total biophysical impacts on LST. When the land is hotter and/or wetter than the atmosphere, a reduced aerodynamic resistance can enhance the sensible and/or latent heat transfer from the land surface to the atmosphere and hence cools the surface (*1*, *28*, *29*). It can also lead to a land surface warming with a reversed heat and water vapor transfer from the atmosphere to the land surface when the land is cooler than the atmosphere, which is reflected in some northern high-latitude regions (fig. S4B). The dominant role of *r*_{a} is robust across results computed from a hybrid of CLM5, MERRA-2, and CMIP5 (fig. S6 and table S2) and is consistent with another independent research using a different model [i.e., MPI-ESM (Max Planck Institute Earth System Model)] (*15*).

Compared to the *r*_{a} pathway, the *r*_{s} pathway plays a secondary role in cooling LST (−0.06 ± 0.09 K m^{2} m^{−2}) (Fig. 2B). In terms of the areal fraction, the *r*_{s} pathway is most significant in 29% of croplands and 19% of grasslands, but least in forests and OWV (8%) (Fig. 2A). This is because although *r*_{s} is more sensitive to changes in LAI than *r*_{a} (fig. S4, H and I), LST is much more sensitive to changes in *r*_{a} (i.e., *r*_{s} does not imply a negligible role of latent heat flux because a change in *r*_{a} can also directly cause a change in latent heat flux (Eq. 4) (*14*–*17*, *30*). Using CLM5’s control and sensitivity runs, we compute the sensitivities of sensible heat flux and latent heat flux to LAI. We find that globally, the sensitivity of latent heat flux (6.9 ± 4.9 W m^{−2} per unit LAI) is larger than that of sensible heat flux (−3.5 ± 4.0 W m^{−2} per unit LAI) (fig. S7). This implies that increases in LAI generally cause the latent heat flux to increase more strongly than the sensible heat flux, which actually shows reductions in most places. The total turbulent flux or available energy (i.e., the sum of sensible and latent heat fluxes) increases with increasing LAI.

The α pathway dominates only in a few places in the Arctic and the Sahel (Fig. 2A). In general, the α pathway leads to a minor warming effect of 0.04 ± 0.07 K m^{2} m^{−2} (Fig. 2B). The finding of albedo-related warming, especially in northern high altitudes, is consistent with previous studies (*2*, *19*, *31*, *32*), but its magnitude is much smaller than the turbulent cooling effects from *r*_{a} and *r*_{s} (Fig. 2). Except for evergreen forests, the increase in LAI reinforces the fraction of absorbed photosynthetically active radiation (fPAR) thus reducing α (*33*), leading to the warming of LST (fig. S4, A and G). However, the increase in LAI in evergreen forests has a negligible α-related impact due to the fact that fPAR is almost saturated in evergreen forests (*34*, *35*) and that the increase in LAI enhances the reflectance in near-infrared radiation, which elevates α (fig. S4, A and G) (*21*, *36*, *37*).

## DISCUSSION

### Causality in the attribution analysis

One highlight of our study is that it treats changes in the LAI as the forcing and changes in the LST as the response, thereby avoiding the potential causality issues involved in the coevolution between LST and LAI. To illustrate this, we use a statistical method following (*10*) to calculate the sensitivity of LST to LAI. The new sensitivity, *10*). However, this new sensitivity is about 10 times larger than the sensitivity diagnosed by the physically based attribution model (

The difference between the two methods is likely caused by the air temperature. In our CLM5 simulations and the TRM method, air temperature is treated as a forcing for the land surface. On the other hand, the regression method does not consider the effects of anomalies of air temperature on anomalies of LST independent from the anomalies of precipitation and incoming shortwave radiation (*10*). As a result, the air temperature effects may be spuriously attributed to the albedo effects in the northern high latitudes given the strong feedbacks in these regions that link the rising temperature with albedo changes. To support this, we analytically computed the sensitivity of LST to air temperature (i.e., _{2} concentration with a minor contribution from the Earth greening (*3*, *38*). We conducted another multiple linear regression with air temperature anomalies as an additional independent variable. The resulting

### Independence of the attributable variables

While the simulated cooling effect from the Earth greening in our study is consistent with another modeling study (*9*), our finding of the dominant role of aerodynamic resistance is in contrast with their conclusion that it plays a negligible role. The main reasons for this difference are their use of 2-m air temperature instead of LST and, equally importantly, the independence assumption made by any attribution methods using first-order Taylor series expansion. By neglecting higher-order and cross-order terms, these attribution methods assume that the attributing factors (such as α, *r*_{a}, and *r*_{s} in our study) are independent of each other (*16*). In (*9*), the authors attributed changes in 2-m air temperature to changes in albedo, aerodynamic resistance, latent heat flux, atmospheric shortwave transmissivity, near-surface air emissivity, and atmospheric circulation. Hence, they implicitly assumed that latent heat flux is independent of aerodynamic resistance. However, the dependence of latent heat flux on aerodynamic resistance is clear from the well-accepted parameterization for latent heat flux (Eq. 4) and has been demonstrated empirically using eddy covariance observations from AmeriFlux (*16*). Their assumption of independence between latent heat flux and aerodynamic resistance may therefore alter the attribution results (*17*).

With regard to the stronger sensitivity of LST to changes in *r*_{a} than in *r*_{s}, this is not inconsistent with a previous theoretical study (*30*) showing that evaporation efficiency is the most important parameter in controlling LST dynamics. First, their estimate of evaporation efficiency is for daytime and relatively moist conditions (*30*), while our study includes both daytime and nighttime and dry conditions at the annual scale. Second, their evaporation efficiency is actually dependent on both aerodynamic resistance and surface resistance due to the use of a different parameterization of latent heat flux (*30*). Our analysis shows that the Earth greening causes large changes in latent heat flux due to the reduced aerodynamic resistance.

## Concluding remarks

In this study, we evaluate the biophysical impacts of Earth greening on LST using an attribution method based on the surface energy balance equation. We find that the widespread Earth greening leads to a cooling effect on LST across the globe at the annual scale, which is predominantly attributed to the decrease in aerodynamic resistance. While these small perturbations in LAI tend to alter turbulent processes more than radiative processes globally, radiative processes remain critical in a small proportion of regions in the Arctic and some sparsely vegetated areas. Last, the TRM method provides a new way to diagnose model outputs and can be further used to evaluate whether these biospheric impacts of vegetation would be amplified or hindered in future climates. If the Earth greening continues, the aerodynamic resistance to turbulent transfer will continue to decrease, resulting in stronger instabilities in the atmospheric boundary layer. In the meantime, surface resistance will also decrease, possibly leading to more water vapor into the atmosphere thus affecting the hydrologic cycle. Whether these effects will be detectable by observations and whether the Earth greening can affect other climatic processes (e.g., extreme events) remain to be investigated.

## MATERIALS AND METHODS

### Model experiments

We conduct offline land model simulations using the CLM5, which is part of the Community Earth System Model Version 2 (CESM2) (*39*), to study the biophysical impacts of Earth greening on the surface climate. In other words, our study object is the land surface, not the coupled land-atmosphere system. Therefore, our CLM5 experiments do not explicitly include the feedbacks of LAI changes on the large-scale climate. Instead, we assume that the ambient atmosphere already considers the impacts of LAI changes and is essentially the forcing of the land surface. The atmospheric forcing is taken from the third phase of the Global Soil Wetness Project [GSWP3 (http://hydro.iis.u-tokyo.ac.jp/GSWP3/)]. All simulations are conducted at 0.47° × 0.63° resolution with a constant atmospheric CO_{2} concentration of 367 parts per million by volume in 2000 to exclude effects from vegetation responses to the rising CO_{2} concentration.

CLM5 explicitly parameterizes the LAI impacts on surface biophysical factors. According to the CLM5 technical note, albedo is influenced by LAI through a two-stream approximation where LAI affects extinction and scattering coefficients (like the Beer-Lambert law) (*39*). For aerodynamic resistance, the displacement height and roughness lengths are both functions of LAI (*22*, *39*). For the surface resistance, CLM5 uses the Medlyn’s model in which surface resistance is affected by LAI through photosynthesis and plant hydraulics (*39*–*41*). Last, emissivity is related to LAI with a negative exponential relationship (*39*).

To obtain the sensitivities of LST to biophysical factors and the sensitivities of biophysical factors to LAI (see Eq. 1), we conduct a suite of simulations for the period of 2000 to 2014 with prescribed LAI obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS):

1. Spin-up run (1990–1999): We prescribe the MODIS LAI climatology, i.e., the LAI values change from month to month but have no interannual variabilities or trends. The spin-up run ensures that the model achieves equilibrium under the current climate.

2. Control run (2000–2014): All configurations are the same as the spin-up run, but the simulation period is from January 2000 to December 2014.

3. Sensitivity runs (2000–2014): Compared to the control run, we perturb the MODIS LAI climatology by ±2%. There are still no interannual variabilities or trends for LAI.

To compare our results to a previous study using multiple linear regression method (

*10*), we conduct a historical LAI run using CLM5:4. Historical LAI run (2000–2013): We prescribe monthly MODIS LAI to CLM5. MODIS LAI is derived from satellite observations, which include interannual variabilities.

### MODIS LAI

The monthly MODIS LAI (0.5° × 0.5°, 2001–2013) at plant functional type (PFT) level is downloaded from the University Corporation for Atmospheric Research (UCAR) data archive for CESM2 (https://svn-ccsm-inputdata.cgd.ucar.edu/trunk/inputdata/lnd/clm2/lai_streams/). For each PFT, we calculate the LAI climatology and the perturbed LAI climatology at the monthly scale.

We use Collection 6 Terra and Aqua MODIS LAI products (MOD15A2H and MYD15A2H, available at https://e4ftl01.cr.usgs.gov) (*33*) to obtain the linear trend in annual average LAI (2000–2014). We check the quality flag for the LAI products to exclude low-quality observations contaminated by clouds, aerosols, shadows, snow, and/or ice. We apply the same preprocess procedure as documented in a previous study (*7*) and resample the dataset to 0.5° × 0.5° to conform to the resolution requirement of CLM5’s inputs. We calculate the trend in annual average LAI using the Mann-Kendall test (R package: https://cran.r-project.org/web/packages/zyp/index.html) at *P* ≤ 0.1. The quality of MODIS LAI was comprehensively evaluated against ground-based measurements and through intercomparisons with other satellite-retrieved LAI products (*42*, *43*).The trend of MODIS LAI is also consistent with AVHRR (Advanced Very-High-Resolution Radiometer) LAI and other LAI datasets proved by a wide range of previous studies (*3*, *7*).

### Diagnosis of the biophysical impacts on LST

As a result of postindustrial human activities, the atmospheric CO_{2} concentration has increased and the land surface has undergone substantial changes. These factors collectively contribute to large-scale climate change and the observed widespread vegetation greening (i.e., LAI change) (*44*). Figure S1 shows that both large-scale climate change (the gray paths) and LAI change (the blue paths) can result in changes in surface biophysical factors (such as albedo, aerodynamic resistance, and surface resistance) (*28*, *45*–*49*), which lastly lead to changes in LST (*14*, *16*), denoted as *46*, *47*). For example, when the air temperature increases, LST would respond even if all the biophysical factors remain the same (*20*). Therefore, only part of the LST change is due to vegetation biophysical impacts

### Sensitivities of LST to biophysical factors

Sensitivities of LST to biophysical factors are analytically computed from the outputs of CLM5 prescribed with LAI climatology using the TRM method. We further examined these sensitivities using reanalysis data (MERRA-2) and outputs from seven fully coupled CMIP5 models (ACCESS1-0, ACCESS1-3, CCSM4, HadGEM2-ES, IPSL-CM5A-LR, IPSL-CM5A-MR, and NorESM1-M). We choose these seven CMIP5 models because only they provide all required inputs for the TRM method (*16*).

We start with the surface energy balance equation, which is expressed as*R*_{n} is the net radiation; *S*_{in} and *L*_{in} are the incoming shortwave and longwave radiation, respectively; α and ε are the albedo and emissivity, respectively; *H* and *LE* are the sensible and latent heat fluxes, respectively; and εσ*T*_{s}^{4} is the outgoing longwave radiation where σ is the Stefan-Boltzmann constant and *T*_{s} is the LST. Further connecting *H* and *LE* with LST through the aerodynamic resistance (*r*_{a}) and surface resistance (*r*_{s}) concepts gives*c*_{p} is the specific heat of air at constant pressure, *L*_{v} is the latent heat of vaporization, *T*_{a} is the air potential temperature, and *q*_{a} is the air-specific humidity. The use of aerodynamic resistance (*r*_{a}) and surface resistance (*r*_{s}) gives rise to the name of the attribution method [TRM method (*14*, *16*, *17*)]. We note that the TRM method calculates the bulk *r*_{a} and *r*_{s} using Eqs. 3 and 4.

Substituting Eqs. 3 and 4 into Eq. 2 yields a nonlinear equation for *T*_{s}, provided that all other variables are given as inputs. This equation is further linearized following previous studies (*14*, *16*, *17*), so that an analytical expression for *T*_{s} can be obtained, as follows*e** is the saturation vapor pressure and *P* is the atmospheric pressure at the surface. With the analytical expression for *T*_{s}, we can take the partial derivatives—

We also calculated the sensitivity of LST to air temperature*16*). We exclude all negative aerodynamic resistance and/or surface resistance (physically meaningless) inferred from Eqs. 3 and 4. We calculate the sensitivities for each year and take the median of these sensitivities from 2000 to 2014 for the consequent analyses.

### Sensitivities of biophysical factors to LAI

We use outputs from the CLM5 control run and sensitivity runs to estimate the sensitivities of biophysical factors to LAI (Eq. 18). For example, the albedo sensitivity to LAI is estimated as*r*_{a}, *r*_{s}, ε, and *G*. As mentioned in the model experiments, there are two sensitivity runs (i.e., perturb the MODIS LAI climatology by ±2%). For each year, we calculate the average sensitivities of biophysical factors to LAI estimated from the two combinations of sensitivity runs and the control run. We take the median sensitivities from 2000 to 2014 for the consequent analyses. Because of a lack of analytical forms of the sensitivities of biophysical factors to LAI, they can only be approximated by offline model simulations to exclude the confounding effects from the changes in atmospheric conditions.

### The robustness of the TRM method and the diagnosed sensitivities

To directly examine the validity of the TRM method, first, we compared the diagnosed

Second, we compare the LST sensitivities to biophysical factors (i.e., *50*), which considers the interactions between lands, atmosphere, and oceans, as well as assimilates observational data (*50*). CMIP5 historical runs are fully coupled simulations of recent past that impose changing conditions consistent with observations, which includes the effects of anthropogenic and volcanic influences and solar activities.

Third, we calculate the

Fourth, we estimate the sensitivity of LST to LAI using the identical multiple linear regression method described in (*10*) (denoted as *T*_{s}, δLAI, δPrecipitation, and δ*SW*_{in} are the annual anomalies from the CLM5 historical LAI run. We note that the MODIS LAI product is the input for CLM5ʼs historical LAI run in satellite phenology mode. Naturally, *b* in Eq. 20. Other coefficients *a*, *c*, and *d* are not discussed here. Besides, we conducted an additional multiple linear regression that includes the interannual variability of air temperature to get a new

### Contribution of each biophysical pathway to ∂ T s bio ∂LAI

For each grid cell, the contribution by each biophysical pathway is defined as the following* _{i}* is the contribution by a particular biophysical pathway, and

### The equivalent energy of ∆ T s bio , LAI

We note that the equivalent energy of *Q*) is a rough estimate that assumes a constant emissivity of unity all over places with significant LAI trends. Thus, *Q* is calculated as*A* is the area of a pixel that varies by latitude; *i* denotes each vegetated pixel with nonzero ^{−8} W m^{−2} K^{−4}); *i*th pixel from 2000 to 2014; *t* is the number of seconds in a year (3.1536×10^{7} s); *N* is the total years of the study period, which is 15; and *n* is the sequential number of year (from 1 to 15).

### MODIS land cover type product

The Collection 5.1 MODIS yearly product provides the land cover information at 0.05° × 0.05° known as MCD12C1 (*51*). The overall accuracy is around 75%; for example, misclassifications are reported between savannas and woody savannas or between cereal croplands and grasslands (*51*–*53*). We aggregate the International Geosphere-Biosphere Programme classification types provided by MCD12C1 into four broad biome types: forests, OWV, grasslands, and croplands. Forests consist of evergreen needleleaf forest, evergreen broadleaf forest, deciduous needleleaf forest, deciduous broadleaf forest, and mixed forest. OWV refers to closed shrublands, open shrublands, and woody savannas. Grasslands include savannas and grasslands. Croplands consist of croplands and croplands/natural vegetation mosaic. There are 12 such global maps, 1 for each year from 2001 to 2012. We refine the 12 maps to 1 map by taking the mode class of each grid cell. Last, we convert the spatial resolution to 0.5° × 0.5° (fig. S2B).

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/47/eabb1981/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**

**Funding:**D.L. and M.H. were supported by the U.S. Department of Energy, Office of Science, as part of research in Multi-Sector Dynamics, Earth and Environmental System Modelling Program. P.G. acknowledges funding from the NSF award number 1734156. R.B.M. acknowledges funding by the Alexander von Humboldt Foundation. We thank Boston University and the National Center for Atmospheric Research (supported primarily by the NSF) for providing supercomputing resources.

**Author contributions:**D.L. and C.C. conceived the methods and experiments. C.C. carried out model simulations and data analyses. C.C. and D.L. wrote the article. All authors contributed comments and critiques on drafts.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All raw data analyzed in this study are publicly available as referenced within the article. CESM2/CLM5 release code is available at www.cesm.ucar.edu/models/cesm2/release_download.html. The trend test package is available at https://cran.r-project.org/web/packages/zyp/index.html. The scripts of the TRM algorithm are available at https://github.com/yaoganchenchi/TRM_Science_Advances. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).