Biophysical impacts of Earth greening largely controlled by aerodynamic resistance.

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 (, ). Satellite records of leaf area index (LAI), a measure of the aboveground vegetation abundance, indicate widespread increasing trends since the 1980s (), due to warming (, ), CO2 fertilization (), and land management (, ). 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 (ra), which characterizes the efficiency of turbulent transfer of heat and water vapor; (iii) surface resistance (rs), 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 (, –). 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 (). 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 (–). 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 (), may suffer causality issues [see, e.g., a recent debate on whether greening is the cause or effect of warming in the boreal zone (–)]. 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 (, , ). It quantifies the sensitivities of LST to Earth greening through biophysical pathways [α, ra, rs, ε, 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 (, , , ), 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 (, , ).

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 (). 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., ). The final expression of the greening-induced LST change is diagnosed as followswhere ∆LAIsat is the satellite-derived long-term trend in LAI (fig. S2A) multiplied by the length of study period (15 years); is the sensitivity of LST to LAI through biophysical pathways; , , , , and are sensitivities of LST to biophysical factors; and , , , , and are sensitivities of biophysical factors to LAI (Materials and Methods).

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 ) indicates a cooling effect due to the Earth greening globally (Fig. 1A)—93% of the global vegetated area shows with an average of −0.36 ± 0.22 K m2 m−2 (mean ± 1 SD, where SD indicates spatial variability). We find that the mean magnitude of is larger in temperate regions (−0.44 K m2 m−2) than those in high-latitude (−0.34 K m2 m−2) and tropical regions (−0.29 K m2 m−2) (Table 1). We further investigated by classifying the global vegetation into forests, other woody vegetation (OWV; such as shrubs and savannas), grasslands, and croplands (fig. S2B). is large in OWV (−0.45 K m2 m−2) and croplands (−0.43 K m2 m−2), followed by grasslands (−0.36 K m2 m−2), but weak in forests (−0.23 K m2 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 increases as precipitation decreases (Table 1). The low in forests is also as expected due to the saturation effect—the same change in LAI would lead to a smaller change in biophysical factors (α, ra, and rs) as well as LST where LAI is high (fig. S4, G to I) (, , ). As a result, the magnitude of in regions with LAI ≤1 m2 m−2 (−0.45 K m2 m−2) is much larger than regions with LAI >4 m2 m−2 (−0.09 K m2 m−2) (Table 1).

Fig. 1

Biophysical effects of the Earth greening on LST.

(A to C) Sensitivity of Ts to LAI, , diagnosed by the TRM method using (A) CLM5 outputs, (B) MERRA-2 and CLM5 outputs, and (C) CMIP5 multi-model ensemble mean (MMEM) and CLM5 outputs. (D) Sensitivity of Ts to LAI, , numerically approximated by CLM5 sensitivity experiments. (E) Sensitivity of Ts to LAI, , estimated by the multiple linear regression method using CLM5 outputs. (F) Changes in Ts due to LAI through biophysical pathways, . All are from 2000 to 2014 except (C) and (E), which are from 2000 to 2005 and 2000 to 2013, respectively.

Table 1

Mean and SD of the biophysical sensitivity of LST to LAI across bioclimatic regimes.

Mean ± 1 SD, where SD indicates the spatial variability. OWV, other woody vegetation.

	Forests	OWV	Grasslands	Croplands	All vegetation
Global	−0.23 ± 0.12	−0.45 ± 0.32	−0.36 ± 0.23	−0.43 ± 0.17	−0.36 ± 0.22
By latitude
>50°S/N	−0.28 ± 0.09	−0.37 ± 0.36	−0.07 ± 0.52	−0.44 ± 0.16	−0.34 ± 0.23
25°S/N–50°S/N	−0.26 ± 0.09	−0.58 ± 0.29	−0.47 ± 0.31	−0.45 ± 0.15	−0.44 ± 0.23
25°S–25°N	−0.11 ± 0.05	−0.45 ± 0.23	−0.33 ± 0.13	−0.37 ± 0.18	−0.29 ± 0.19
By LAI
LAI < 1 m2 m−2	−0.38 ± 0.15	−0.49 ± 0.38	−0.37 ± 0.37	−0.46 ± 0.17	−0.45 ± 0.31
LAI ∈1–4 m2 m−2	−0.25 ± 0.08	−0.33 ± 0.11	−0.33 ± 0.11	−0.39 ± 0.14	−0.30 ± 0.11
LAI > 4 m2 m−2	−0.09 ± 0.02	−0.13 ± 0.03	−0.12 ± 0.02	−0.11 ± 0.02	−0.09 ± 0.02
By annual total precipitation (ATP)
ATP < 900 mm	−0.29 ± 0.09	−0.47 ± 0.36	−0.35 ± 0.30	−0.43 ± 0.17	−0.40 ± 0.26
ATP ∈ 900–2000 mm	−0.21 ± 0.10	−0.35 ± 0.13	−0.37 ± 0.15	−0.43 ± 0.16	−0.33 ± 0.16
ATP > 2000 mm	−0.10 ± 0.04	−0.31 ± 0.17	−0.31 ± 0.19	−0.33 ± 0.18	−0.12 ± 0.06

Multiple lines of evidence show that our diagnosed from CLM5 outputs using the TRM method is robust. First, (Fig. 1A) is almost identical to its counterpart (i.e., ; Materials and Methods) approximated using CLM5 sensitivity experiments by perturbing LAI climatology (Fig. 1D). This directly demonstrates that the TRM method can successfully capture the nonlinear response of LST to LAI simulated by CLM5. Second, given that there have been studies reporting CLM’s biases in modeling sensible and latent heat fluxes (, ), 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., ), and we find that the results agree well with those diagnosed from CLM5 outputs (Fig. 1, B and C, and fig. S5). Third, we also use these LST sensitivities to biophysical factors computed from MERRA-2 and CMIP5 to calculate (fig. S6). The coefficient of variation (CV) of these is 15%, indicating good agreement among models. We also break down into three individual biophysical pathways. Values of CV for each pathway are 44% (albedo), 10% (aerodynamic resistance), and 22% (surface resistance). The relatively high CV in the albedo pathway turns out to be unimportant because the absolute contribution from the albedo pathway to is much smaller than those from aerodynamic and surface resistances (fig. S6). We stress that different inputs for TRM agree in terms of their signs for all biophysical pathways, including the albedo pathway (fig. S6).

Earth greening cools LST

The change of LST due to Earth greening (i.e., ) is the product of sensitivity and the long-term change in LAI (i.e., ∆LAIsat) (Eq. 1). Our results show that is −0.056 ± 0.046 K over regions with statistically significant LAI trends, and the area of LAI-induced cooling (30%) is six times of that of LAI-induced warming (5%) (Fig. 1F). The magnitude of estimated here is comparable to that reported by a previous study using coupled land-atmosphere simulations (). 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 from 2000 to 2014 to changes in energy (2.97 × 1021 J), which is more than five times the world’s total primary energy supply in 2015 (5.71 × 1020 J) (). 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 × 1020 J) (), given that one-third of the global vegetated area shows statistically significant trends in LAI ().

Because the change in LST is the product of the sensitivity, , and the change in LAI, ∆LAIsat (Eq. 1), high sensitivity of LST to LAI does not necessarily imply a large change in LST, and vice versa. For example, is large in semiarid regions and/or highlands, such as in the western United States, western China, Australia, South Africa, and Argentina (Fig. 1A). However, these regions have small ∆LAIsat so that is suppressed (Fig. 1F and fig. S2A). A similar situation exists in the Arctic where is positive and large (Fig. 1A), yet is negligible because ∆LAIsat is small (Fig. 1F and fig. S2A). In contrast, in subtropical to temperate regions of the Northern Hemisphere (eastern China, India, and Eastern Europe), is large (Fig. 1F) due to the strong (Table 1) and ∆LAIsat (fig. S2A).

The dominant role of aerodynamic resistance

We quantify the contributions from different biophysical factors to , including α, ra, rs, ε, and G (Eq. 1). At the annual scale, ra 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 (, ). In CLM5, the ra pathway dominates 82% of the global vegetated area, followed by the rs 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 ra pathway plays a more prominent role than α and rs, while the dominant fraction of ra 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 m2 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 (, , ). 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 ra 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)] ().

Fig. 2

Dominant surface biophysical factors in regulating at the annual scale diagnosed from CLM5 outputs.

(A) Map of dominant factors for . Orange, yellow, and green represent the dominance of α, ra, and rs, respectively. The inset shows the areal fraction of dominant factors by biome type. FO, forests; OWV, other woody vegetation; GR, grasslands; CR, croplands; All, all vegetation. (B) Attribution of with surface biophysical factors. Results are presented in boxplot, and the additional diamonds indicate the mean.

Compared to the ra pathway, the rs pathway plays a secondary role in cooling LST (−0.06 ± 0.09 K m2 m−2) (Fig. 2B). In terms of the areal fraction, the rs 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 rs is more sensitive to changes in LAI than ra (fig. S4, H and I), LST is much more sensitive to changes in ra (i.e., ) in CLM5 (fig. S4, E and F) as well as in MERRA-2 and CMIP5 (fig. S5, K and L). We stress that a minor role of rs does not imply a negligible role of latent heat flux because a change in ra can also directly cause a change in latent heat flux (Eq. 4) (–, ). 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 m2 m−2 (Fig. 2B). The finding of albedo-related warming, especially in northern high altitudes, is consistent with previous studies (, , , ), but its magnitude is much smaller than the turbulent cooling effects from ra and rs (Fig. 2). Except for evergreen forests, the increase in LAI reinforces the fraction of absorbed photosynthetically active radiation (fPAR) thus reducing α (), 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 (, ) and that the increase in LAI enhances the reflectance in near-infrared radiation, which elevates α (fig. S4, A and G) (, , ).

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 () to calculate the sensitivity of LST to LAI. The new sensitivity, , is calculated using multiple linear regression of LST anomalies against LAI, precipitation, and incoming shortwave radiation anomalies from CLM5 (Fig. 1E). This new sensitivity can almost reproduce the results of (). However, this new sensitivity is about 10 times larger than the sensitivity diagnosed by the physically based attribution model (), and the two have opposite signs in the northern high latitudes (Fig. 1, A and E). Since both and are computed from CLM5 outputs, they are expected to agree with that evaluated from CLM5 sensitivity runs with prescribed LAI perturbations (; Fig. 1D and table S1). However, only the former is consistent with this expectation.

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 (). 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., ) using the TRM method, which shows strong positive values over the northern high-latitude regions (fig. S8A). This implies that the LST and air temperature in these regions are strongly and positively correlated, and the rising air temperature inevitably leads to increased LST. The rising air temperature in these regions is more likely caused by the increasing atmospheric CO2 concentration with a minor contribution from the Earth greening (, ). We conducted another multiple linear regression with air temperature anomalies as an additional independent variable. The resulting becomes negative in most places albeit with a different magnitude compared to the TRM method (fig. S8B), which suggests that results from the regression method are highly dependent on the selected “independent” forcing variables.

Independence of the attributable variables

While the simulated cooling effect from the Earth greening in our study is consistent with another modeling study (), 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 α, ra, and rs in our study) are independent of each other (). In (), 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 (). Their assumption of independence between latent heat flux and aerodynamic resistance may therefore alter the attribution results ().

With regard to the stronger sensitivity of LST to changes in ra than in rs, this is not inconsistent with a previous theoretical study () 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 (), 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 (). 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) (), 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 CO2 concentration of 367 parts per million by volume in 2000 to exclude effects from vegetation responses to the rising CO2 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) (). For aerodynamic resistance, the displacement height and roughness lengths are both functions of LAI (, ). For the surface resistance, CLM5 uses the Medlyn’s model in which surface resistance is affected by LAI through photosynthesis and plant hydraulics (–). Last, emissivity is related to LAI with a negative exponential relationship ().

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):

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.

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.

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 (), we conduct a historical LAI run using CLM5:

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) () 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 () 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 (, ).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 (, ).

Diagnosis of the biophysical impacts on LST

As a result of postindustrial human activities, the atmospheric CO2 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) (). 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) (, –), which lastly lead to changes in LST (, ), denoted as (the gray and red paths) and (the blue and red paths), respectively. In addition, large-scale climate change can cause changes in LST without altering the biophysical factors (the pink path, denoted as ; fig. S1) (, ). For example, when the air temperature increases, LST would respond even if all the biophysical factors remain the same (). Therefore, only part of the LST change is due to vegetation biophysical impacts . The full expression of is Eq. 1 in the main text.

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

We start with the surface energy balance equation, which is expressed aswhere Rn is the net radiation; Sin and Lin 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 εσTs4 is the outgoing longwave radiation where σ is the Stefan-Boltzmann constant and Ts is the LST. Further connecting H and LE with LST through the aerodynamic resistance (ra) and surface resistance (rs) concepts giveswhere ρ is the air density, cp is the specific heat of air at constant pressure, Lv is the latent heat of vaporization, Ta is the air potential temperature, and qa is the air-specific humidity. The use of aerodynamic resistance (ra) and surface resistance (rs) gives rise to the name of the attribution method [TRM method (, , )]. We note that the TRM method calculates the bulk ra and rs using Eqs. 3 and 4.

Substituting Eqs. 3 and 4 into Eq. 2 yields a nonlinear equation for Ts, provided that all other variables are given as inputs. This equation is further linearized following previous studies (, , ), so that an analytical expression for Ts can be obtained, as followswhereand e* is the saturation vapor pressure and P is the atmospheric pressure at the surface. With the analytical expression for Ts, we can take the partial derivatives—, , ,, and to obtain sensitivities of LST to biophysical factors, as follows

We also calculated the sensitivity of LST to air temperaturewhere , , , and . To use the TRM attribution method, the required input variables are from CLM5 control run, which include incoming shortwave radiation, outgoing shortwave radiation, incoming longwave radiation, sensible heat flux, latent heat flux, ground heat flux, emissivity, surface pressure, and air temperature, specific humidity, and air density at the lowest atmospheric model level (about 30 m above the land surface) (). 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 aswhere the subscripts “sen” and “ctl” indicate sensitivity and control simulations, respectively. Note that α can be replaced with ra, rs, ε, 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 (by TRM) against the sensitivity () directly calculated from CLM5 control and sensitivity runs. Similar to the biophysical sensitivities to LAI, the can be calculated as

Second, we compare the LST sensitivities to biophysical factors (i.e., , , and ) calculated from CLM5 to those (i) from MERRA-2 (0.625° × 0.5°, 2000–2017) and (ii) from seven CMIP5 historical runs (resolution varies by model, 2000–2005). This can address whether the analytically analyzed sensitivities (i.e., , , and ) diagnosed from offline CLM5 simulations are consistent with those diagnosed from reanalysis data and fully coupled CMIP5 simulations. MERRA-2 is the latest version of reanalysis data produced by the NASA’s Global Modelling and Assimilation Office (), which considers the interactions between lands, atmosphere, and oceans, as well as assimilates observational data (). 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 using a combination of MERRA-2, CMIP5, and CLM5. That is, , , and are estimated from MERRA-2 and CMIP5, and , , and are estimated from CLM5. In other words, we compared the hybrid to the estimated purely from CLM5 outputs.

Fourth, we estimate the sensitivity of LST to LAI using the identical multiple linear regression method described in () (denoted as )where δTs, δLAI, δPrecipitation, and δSWin 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, is the slope 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 to compare with the estimated by Eq. 20.

Contribution of each biophysical pathway to

For each grid cell, the contribution by each biophysical pathway is defined as the followingwhere Contribution is the contribution by a particular biophysical pathway, and .

The equivalent energy of

We note that the equivalent energy of (denoted as Q) is a rough estimate that assumes a constant emissivity of unity all over places with significant LAI trends. Thus, Q is calculated aswhere A is the area of a pixel that varies by latitude; i denotes each vegetated pixel with nonzero ; σ is the Stefan-Boltzmann constant (5.67×10−8 W m−2 K−4); is the climatology LST of the ith pixel from 2000 to 2014; t is the number of seconds in a year (3.1536×107 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 (). The overall accuracy is around 75%; for example, misclassifications are reported between savannas and woody savannas or between cereal croplands and grasslands (–). 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).