# Biophysical impacts of northern vegetation changes on seasonal warming patterns

Xu Lian, Sujong Jeong, Chang-Eui Park, Hao Xu, Laurent Z. X. Li, Tao

Wang, Pierre Gentine, Josep Peñuelas, Shilong Piao

► To cite this version:

Xu Lian, Sujong Jeong, Chang-Eui Park, Hao Xu, Laurent Z. X. Li, et al. Biophysical impacts of northern vegetation changes on seasonal warming patterns. Nature Communications, 2022, 13, 10.1038/s41467-022-31671-z. insu-03847110

HAL Id: insu-03847110

https://insu.hal.science/insu-03847110

Submitted on 11 Nov 2022

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>

ARTICLE

https://doi.org/10.1038/s41467-022-31671-z

OPEN

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

Check for updates

# Biophysical impacts of northern vegetation changes on seasonal warming patterns

Xu Lian<sup>1,2</sup>, Sujong Jeong<sup>3,4✉</sup>, Chang-Eui Park<sup>5</sup>, Hao Xu<sup>1</sup>, Laurent Z. X. Li<sup>6</sup>, Tao Wang<sup>7</sup>, Pierre Gentine<sup>2,8</sup>, Josep Peñuelas<sup>9,10</sup> & Shilong Piao<sup>1,7✉</sup>

The seasonal greening of Northern Hemisphere (NH) ecosystems, due to extended growing periods and enhanced photosynthetic activity, could modify near-surface warming by perturbing land-atmosphere energy exchanges, yet this biophysical control on warming seasonality is underexplored. By performing experiments with a coupled land-atmosphere model, here we show that summer greening effectively dampens NH warming by -0.15 ± 0.03 °C for 1982-2014 due to enhanced evapotranspiration. However, greening generates weak temperature changes in spring (+0.02 ± 0.06 °C) and autumn (-0.05 ± 0.05 °C), because the evaporative cooling is counterbalanced by radiative warming from albedo and water vapor feedbacks. The dwindling evaporative cooling towards cool seasons is also supported by state-of-the-art Earth system models. Moreover, greening-triggered energy imbalance is propagated forward by atmospheric circulation to subsequent seasons and causes sizable time-lagged climate effects. Overall, greening makes winter warmer and summer cooler, attenuating the seasonal amplitude of NH temperature. These findings demonstrate complex tradeoffs and linkages of vegetation-climate feedbacks among seasons.

<footnote n="1">¹Sino-French Institute for Earth System Science, College of Urban and Environmental Sciences, Peking University, Beijing, China. ²Department of Earth and Environmental Engineering, Columbia University, New York, NY, USA. ³Department of Environmental Planning, Graduate School of Environmental Studies, Seoul National University, Seoul, Republic of Korea. ⁴Environmental Planning Institute, Seoul National University, Seoul, Republic of Korea. ⁵Center for Sustainable Environment Research, Korea Institute of Science and Technology, Seoul, Republic of Korea. ⁶Laboratoire de Météorologie Dynamique, CNRS, Sorbonne Université, Ecole Normale Supérieure, Ecole Polytechnique, Paris, France. ⁷State Key Laboratory of Tibetan Plateau Earth System, Resources and Environment (TPESRE), Institute of Tibetan Plateau Research, Chinese Academy of Sciences, Beijing, China. ⁸Earth Institute, Columbia University, New York, NY, USA. ⁹CREAF, Cerdanyola del Valles, Barcelona, Catalonia, Spain. ¹⁰CSIC, Global Ecology Unit CREAF-CSIC-UAB, Bellaterra, Barcelona, Catalonia, Spain.
✉email: sujong@snu.ac.k; slpiao@pku.edu.cn</footnote>
As air temperature rises, the phenological cycle of Northern Hemisphere (NH) ecosystems is shifting progressively towards earlier leaf emergence and delayed leaf senescence, which leads to rapid lengthening of the active growing season¹⁻⁴. Remote sensing techniques have also clearly documented a widespread increase in leaf area and enhanced vegetation activity throughout the growing season, a phenomenon known as “vegetation greening”⁵⁻⁸. In turn, the altered phenology and structure of terrestrial plants could modulate the warming rate by changing the seasonal cycles of carbon, energy, and water at the land-atmosphere interface⁴,⁹. However, current knowledge about this climatic effect is limited mostly to enhanced carbon sequestration caused by a lengthened growing period and accelerated photosynthesis rates (biochemical feedbacks)¹,⁴,¹⁰. The role of vegetation changes in the seasonal budget of surface energy fluxes (biophysical feedbacks)¹¹,¹² has been studied far less extensively, despite their strong capacity to affect regional to global warming over annual or longer timescales¹³⁻¹⁵.

Vegetation biophysics have long been recognized as a key regulator of seasonal air temperature climatology. For example, in early spring, the rate of air temperature increase rapidly decreases after leaf unfolding (typically for deciduous forests) due to increased transpiration after leaf-out that can effectively cool the leaf surface¹²,¹⁶,¹⁷. The critical role of plants in local temperature seasonality suggests that greening will alter the seasonality of NH warming at annual to decadal timescales. Vegetation greening also affects climate by interacting with other land-surface (e.g., snow or soil moisture) and atmospheric (e.g., water vapor, cloud, and circulation) processes¹¹,¹³,¹⁸, the effects of which vary both geographically and seasonally⁴. For example, during snow-affected seasons/areas, vegetation greening markedly decreases snow cover and surface albedo due to a strong masking effect of darker canopies¹¹,¹⁸⁻²⁰, and the surface darkening would warm the surface that further accelerates snow losses²¹,²². In addition to affecting climate locally and transiently, vegetation biophysical processes summarized above can potentially trigger climate anomalies beyond the greening region (non-local feedbacks) and season (inter-seasonal feedbacks)²³⁻²⁵, by redistributing available energy over space and time. In particular, potential time-lagged temperature responses are hidden in those annually aggregated climatic consequences of greening¹³, and to some degree affect the seasonally asymmetric characteristics of NH warming²⁶⁻²⁸. However, our understanding of such effects and underlying processes remains rudimentary so far.

Observation-driven assessments compare the temperature difference between undisturbed forests and their neighboring non-forest counterparts as a proxy for temperature responses to forest gain, assuming the same climate background¹⁴,²⁹,³⁰. This space-for-time substitution is effective for quantifying localized and transient temperature feedbacks¹⁴,²⁹,³⁰, yet overlooks indirect atmospheric feedbacks due to, for example, water vapor, cloud, and circulation changes³¹,³². Hence, it is intrinsically flawed to detect non-local and/or inter-seasonal signals pertained to seasonal vegetation-climate feedbacks. Coupled land-atmosphere global climate models (GCM) offer a useful tool to detect causality in complex systems, e.g., the vegetation-climate interactions, though with uncertainties from a simplified representation of real-world processes³³. GCMs link vegetation dynamics to climate changes by parameterizing the leaf area index (LAI) that determines land-surface boundary conditions and controls vegetation-atmospheric exchanges of water and energy¹³,³⁴, and explicitly represent large-scale changes of atmospheric physical processes.

In this study, we use a series of simulations with the IPSL-CM (Institute Pierre Simon Laplace coupled model, see Methods) GCM³⁵ to examine the biophysical effects of the observed widespread NH (defined here as 25–90°N) greening on seasonal

warming patterns. Our key finding is that the cooling effect of NH greening is most prominent in summer, which, however, diminishes or even turns to net warming in cool seasons due to weakened evaporative capacity and enhanced albedo and water vapor feedbacks. We also show sizable time-lagged climate effects that greening in one season affects the temperature in the following ones by perturbing atmospheric circulation patterns. Such biophysical control of greening ultimately leads to a deviation of the seasonal warming trajectory from that originally driven by greenhouse warming.

## Results and discussion

**Coupled model experiments for detecting vegetation-climate feedback.** We quantified changes of near-surface (2-m) air temperature ($T_a$) in response to the observed NH greening for all active growing seasons during 1982–2014 using IPSL-CM. We defined the three growing seasons (spring, summer, and autumn) across the entire NH domain as periods of March-April-May (MAM), June-July-August (JJA), and September-October-November (SON), respectively. For each season, a pair of transient numerical experiments was performed by modifying LAI: a dynamic vegetation experiment (SCE) forced by annually and seasonally varying LAI from satellite observations³⁶, and three seasonal control experiments ($LAI_{CTL}^{MAM}$, $LAI_{CTL}^{JJA}$, and $LAI_{CTL}^{SON}$ for MAM, JJA, and SON, respectively) forced by annually varying LAI for all seasons, except in the season of interest when the LAI was fixed to the climatological conditions observed during 1982–2014 (Fig. S1). For all experiments, other boundary conditions, including sea surface temperature (SST), sea ice fraction (SIC), and atmospheric $CO_2$ concentrations, were kept consistent (Methods). Therefore, differences between SCE and the control experiments characterized the effects of the observed LAI changes on $T_a$ (hereafter denoted as $\Delta T_a$), both intra- and inter-seasonally. Multimember paired ensembles were generated for each coupled model experiment by performing 30 repeated runs but with different initial conditions (see Methods).

The capacity of the IPSL-CM GCM for simulating the seasonal variations and spatial patterns of $T_a$ was assessed by comparing the SCE simulation results with the observation-based $T_a$ data (Methods). Throughout most of the growing season (May to October), the SCE simulation well reproduced the increasing trend and interannual variability of the NH land mean $T_a$ observed during 1982–2014 (Fig. S2). Observational data showed that the strongest NH warming occurred in early spring (March and April) and late autumn (November). However, the SCE simulation failed to capture the exceptionally strong warming during the transitional seasons, leading to the underestimation of the annual mean warming trend (SCE: $0.237 \pm 0.024$ °C decade⁻¹; observed: $0.362 \pm 0.048$ °C decade⁻¹). This underestimation stemmed from a negative bias in the increase of downwelling shortwave radiation, possibly due to an absence of short-lived forcing and bias in the cloud systems³⁷. Overall, the SCE reproduced the geographical patterns of seasonal warming reasonably well (Fig. S3), which strengthened our confidence in the model projections. Notably, it successfully captured the observed amplified warming over pan-arctic and semi-arid regions, as well as the few cases of regional cooling, such as that over northwestern North America during MAM (Fig. S3).

**Intra-seasonal temperature responses to NH LAI changes.** For the period from 1982 to 2014, satellite-retrieved LAI showed statistically significant increasing trends ($p < 0.05$) during all growing seasons, with the most substantial increase observed in JJA (0.046 ± 0.007 $m^2 m^{-2}$ decade⁻¹), followed by SON (0.027 ± 0.005 $m^2 m^{-2}$ decade⁻¹) and MAM (0.019 ± 0.005 $m^2 m^{-2}$ decade⁻¹)⁵ (Fig. 1a).
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Intra- and inter-seasonal temperature responses to leaf area index (LAI) changes. **a** Monthly trends (shadings) of Northern Hemisphere (NH) mean LAI during 1982–2014 used as input to the seasonal simulations. The dashed curve and transparent bars indicate trends of monthly LAI and seasonally aggregated LAI values, respectively. **b** Linear trends of $T_a$ driven by LAI changes within the same season (intra-seasonal) and other growing seasons (inter-seasonal). Error bars in **a**, **b** indicate uncertainty ranges [1 – standard deviation (SD)]. **c** Monthly trends of LAI-induced air temperature changes ($\Delta T_a$), with red and blue shadings representing positive and negative trends, respectively. The bottom panel shows the overall $\Delta T_a$ trends induced by LAI changes in all growing seasons, calculated as the sum of $\Delta T_a$ trends from the three seasonal runs shown separately in the above panels. ***$p < 0.01$; **$p < 0.05$; *$p < 0.1$; n.s., $p > 0.1$. MAM March-April-May, JJA June-July-August, SON September-October-November, DJF December-January-February.</figcaption></figure>

Although forced identically by observed greening patterns, our GCM estimated highly divergent $T_a$ responses across seasons, showing non-significant MAM warming ($0.005 \pm 0.018$ °C decade$^{-1}$, $p > 0.1$), strong and significant JJA cooling ($-0.044 \pm 0.008$ °C decade$^{-1}$, $p < 0.05$), and non-significant SON cooling ($-0.014 \pm 0.016$ °C decade$^{-1}$, $p > 0.1$) (intra-seasonal feedbacks shown in Fig. 1b). The LAI-induced JJA $T_a$ trend was equivalent to cooling of $-0.15 \pm 0.03$ °C in JJA over the study period, offsetting the overall SCE-simulated near-surface air warming over this period by ~12.5%. This strong JJA cooling was further supported by a significant negative correlation ($r = -0.64$, $p < 0.05$) between JJA LAI and its induced $T_a$ anomalies (Fig. S4b), with greener summers (higher JJA LAI) co-occurring with more negative $\Delta T_a$. However, no such correlation was statistically detectable in MAM ($r = -0.14$, $p > 0.1$) or SON ($r = 0.07$, $p > 0.1$) (Fig. S4a, c), during which the LAI-induced changes accounted for only 1.3% (MAM) and ~3.2% (SON) of the concurrent greenhouse warming. We also verified the robustness of our results by performing equilibrium experiments with an independent model, the NCAR Community Atmosphere Model coupled with Community Land Model (CAM-CLM, Methods). Indeed, this model generated a similarly strong LAI-induced cooling in JJA ($-0.18$ °C, $p < 0.05$ based on a two-sample *t*-test), which shifted to be much weaker in MAM ($-0.02$ °C, $p > 0.1$) and SON ($-0.05$ °C, $p > 0.1$) (Fig. S5).

Figure 2 shows the spatial distribution of the intra-seasonal $T_a$ responses to LAI changes based on IPSL-CM. In the JJA experiment (SCE – LAI<sup>JJA</sup><sub>CTL</sub>), northern areas where $\Delta T_a$ decreased the most corresponded well to those where LAI showed the strongest greening, such as Europe, Russia, and the central U.S. (Fig. 2b, e). Similarly, increasing $\Delta T_a$ occurred in land areas that experienced the strongest browning, such as the southern U.S., Alaska, and semi-arid Asia

(Fig. 2b, e). However, this spatial consistency of greening-cooling (or browning-warming) was not detected in most northern areas during MAM or SON (Fig. 2a, c, d, f). During these two seasons, though widespread areas showed significant greening (MAM: 42.0%; SON: 40.3%), only a small fraction (<5%) of northern lands exhibited significant trends of extra cooling or warming (Fig. 2a, c, d, f). In regions of Europe, East Asia, and the eastern U.S., the regional surface cooling in MAM and SON was also insignificant and discernably weaker than that in JJA, albeit with comparably strong greening signals across the three seasons. CAM-CLM also broadly agreed with the spatial patterns of signs despite varying magnitudes (Fig. S6). On one hand, these results imply a lower capacity for spring and autumn vegetation changes to influence near-surface climate, with potentially feedback processes different from summer periods. On the other hand, process-based uncertainties are more substantial in these relatively cold seasons associated with a model deficiency in representing vegetation–snow–moisture interactions.

**Inter-seasonal temperature responses to NH LAI changes.** The IPSL-CM model also produced particularly strong time-lagged $T_a$ responses to greening in previous growing seasons. For example, MAM greening-induced cooling in late spring carried over to early summer; JJA greening-induced cooling extended throughout the entire growing season, although cooling in autumn and spring of the next year was not as strong as that in summer; and SON LAI changes caused anomalously strong cooling in the following spring (Fig. 1c). In all growing seasons, NH greening consistently resulted in net warming during winter (December to February, DJF). For each focused season, the overall time-lagged response of $T_a$ to previous-season greening was calculated by summing up the
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Spatial patterns of intra-seasonal temperature responses to leaf area index (LAI) changes. Maps of linear trends in seasonal LAI (a-c) and the induced air temperature changes (ΔT_a) in the corresponding season (d-f) during 1982-2014. ΔT_a values reflect the differences in results between the dynamic vegetation experiment (SCE) and the three seasonal control experiments (LAI<sup>MAM</sup><sub>CTL</sub>, LAI<sup>JJA</sup><sub>CTL</sub>, and LAI<sup>SON</sup><sub>CTL</sub>). Stippling indicates areas where the trend is statistically significant (*p* < 0.05). Pixels with climatological LAI < 0.1 are masked. MAM March-April-May, JJA June-July-August, SON September-October-November, DJF December-January-February.</figcaption></figure>

GCM-simulated T_a changes in response to greening in the two other growing seasons. At the scale of the NH, the GCM estimated ΔT_a trends of -0.057 ± 0.029 °C decade⁻¹ (*p* = 0.06; MAM), -0.016 ± 0.013 °C decade⁻¹ (*p* = 0.21; JJA), and -0.029 ± 0.019 °C decade⁻¹ (*p* = 0.13; SON) in response to previous-season greening (Fig. 1b), corresponding to cooling of 0.19 ± 0.10 °C, 0.05 ± 0.04 °C, and 0.10 ± 0.06 °C, respectively during 1982–2014. In JJA, this cooling signal caused by previous greening amplified the transient cooling caused by JJA greening by ~36%. In MAM and SON, this inter-seasonal signal exceeded that caused by LAI changes within the same season (Fig. 1b), as also confirmed by the CAM-CLM results (Fig. S5), thereby dominating their overall T_a responses to LAI changes for the entire growing season. In both models, however, we did not detect the inter-seasonal signal with a very high statistical confidence level, due to trade-offs among regions and processes driving this carry-over effect.

Examination of the intra- and inter-seasonal biophysical
feedbacks together revealed that growing season vegetation
greening had a discernable effect on the seasonal trajectory of
the NH land mean T_a (Fig. 1c, bottom panel). Seasonal greening
overall contributed to net warming in winter and net cooling
throughout the entire growing season. By simultaneously causing
warmer winters and cooler summers, NH greening reduced the
amplitude of the T_a annual cycle (SAT, defined as July T_a minus
January T_a) at an average rate of -0.14 ± 0.07 °C decade⁻¹
(*p* < 0.1) (Fig. 3), signaling to weaken T_a seasonality. Several
previous studies reported that land SAT decreased throughout the
twentieth century<sup>26-28</sup>. However, during the more recent
1982-2014 period, real-world observations and SCE results
consistently showed a trend of leveling-off or a non-significant

<figure><img src="image_6.png" /><figcaption>**Fig. 3** Vegetation-induced changes in the annual temperature range.</figcaption></figure>

Anomalies of observed (brown, from the Princeton Global Meteorological
Forcing product) and IPSL-CM-simulated (purple, from the dynamic
vegetation experiment SCE) annual ranges of Northern Hemisphere (NH)
land mean air temperature (T_a). The bold brown curve indicates the 10-year
running mean of the observed annual T_a range. The green curve indicates
the simulated annual T_a range induced by leaf area index (LAI) changes in
all growing seasons combined. \* *p* < 0.1; n.s., *p* > 0.1.

increase (+0.11 °C decade⁻¹) (Fig. 3), indicating that greening
diminished the recent increase in the seasonal amplitude of T_a in
a warmer climate. In other words, in the absence of vegetation-
climate feedbacks, northern lands would have experienced a
<figure><img src="image_7.png" /><figcaption>**Fig. 4 Biophysical processes underlying intra-seasonal temperature responses to greening. a–c** Linear trends in surface radiative forcing associated with changes in latent heat flux (λE), surface albedo (α), aerodynamic resistance (r_a), shortwave radiation (S_dn), and air emissivity (ε_a) induced by leaf area index (LAI) changes within the same season. **d–f** Latitudinal variation in surface radiative forcing associated with different biophysical processes. “All” indicates the sum of all surface radiative forcings; “R” and “NR” indicate surface radiative forcing associated with radiative and non-radiative processes, respectively. ***p < 0.01; **p < 0.05; *p < 0.1; n.s., p > 0.1.** MAM March-April-May, JJA June-July-August, SON September-October-November, DJF December-January-February.</figcaption></figure>

much faster SAT increase, by about +0.25 °C decade⁻¹, over the past 33 years.

Driving processes of the seasonal vegetation-climate feedback.

An intriguing question is why the detected vegetation biophysical
feedbacks on T_a vary in sign and magnitude among seasons. This
variation is partly explained by the differing greening trends
among seasons, with the most significant greening and cooling
co-occurring in summer (Fig. 1a, b). Moreover, the major bio-
physical processes that govern the surface energy budget may also
diverge seasonally. To fully understand the mechanisms of sea-
sonal dependence among LAI-T_a feedbacks, we decomposed the
modeled T_a response to LAI changes into contributions from
different surface physical properties/processes, including surface
albedo (α), evapotranspiration (λE or ET), surface incoming
shortwave radiation (S_dn), air longwave radiative emissivity (ε_a),
and aerodynamic resistance (r_a) (see Methods). Together, these
processes well reproduced the modeled LAI-induced net change
in surface radiative forcing related to changes in T_a (Fig. 4a–c).

During all growing seasons, the ET-related change in surface
radiative forcing was the greatest among the surface processes
considered (Fig. 4a–c). The ET enhancement was paralleled with
a proportional decrease in sensible heat (Fig. S7a–c). The less
available surface energy dissipated as latent heat (Fig. S7d) favors
the evaporative cooling of the surface. Geographically, regions
with ET enhancement (Fig. 4d–f and S8b, h, n) also corresponded
well to those of satellite-observed greening (Fig. 2a–c), confirming
localized vegetation-temperature feedbacks through this key
process¹³,²⁹. Among the three growing seasons, greening
generated a greater increase in the fraction of net solar radiation
partitioned into ET (Fig. S7d), which resulted in a greater cooling

in JJA (−0.64 ± 0.09 W m⁻² decade⁻¹, *p* < 0.05) than in MAM
(−0.23 ± 0.05 W m⁻² decade⁻¹, *p* < 0.05) and SON
(−0.15 ± 0.03 W m⁻² decade⁻¹, *p* < 0.05) (Fig. 4a–c). In JJA,
plants often consume a greater amount of soil water to keep up
with the greater water demand. For example, the ratio of
transpiration to total ecosystem ET (T/ET)—a measure of
vegetation control on surface climate through evaporative
cooling³⁸—was ~50% greater in JJA than in MAM/SON
(Fig. S9a), which resulted in larger positive sensitivity of land
ET to LAI translated to stronger surface cooling in JJA (Fig. S9a).
The stronger evaporative cooling capacity of JJA greening than
other growing seasons was also supported by 33 state-of-the-art
Earth system models (Fig. S9b), suggesting that the seasonally
varying LAI-T_a feedbacks are not dependent on our specific
model use. Interestingly, within both MAM and SON, the LAI-T_a
feedbacks shifted progressively from warming-dominated in
colder months (i.e., March and November) to cooling-
dominated in warmer months (i.e., May and September) (Fig. 1c).
This phenomenon can be partly explained by the
land-atmosphere theory that latent heat of vaporization is a
more efficient cooling mechanism at higher temperatures³⁹.

In addition to affecting how the land surface dissipates absorbed
energy, greening also affects how much energy the surface absorbs
(i.e., due to changes in α, S_dn, and ε_a). These radiative processes
played a secondary, but often non-negligible, role, together
increasing the surface radiative forcing at an average rate of
0.25 W m⁻² decade⁻¹ (*p* < 0.05; MAM), 0.18 W m⁻² decade⁻¹
(*p* < 0.05; JJA), and 0.05 W m⁻² decade⁻¹ (*p* < 0.05; SON)
(Fig. 4a–c). During MAM when ET is constrained by relatively
low temperatures, this radiative warming (as the combined effect
of α, S_dn, and ε_a) overrode evaporative cooling and caused net
surface warming (Fig. 4a). Increased ε_a contributed the most to
this warming (Fig. 4a, d), as greening increased atmospheric water vapor content, which warmed the air by trapping extra longwave radiation<sup>13,21</sup>. Decreased α also significantly increased surface radiative forcing, particularly over greening and boreal/arctic regions, and throughout the seasonal cycle, such albedo-related warming peaked during MAM in the presence of snow cover<sup>11,18</sup> (Fig. 4a, d and S8c). In JJA, radiative warming, as a result of decreased α and increased ε<sub>a</sub>, offset only 26% of evaporative cooling, causing net surface cooling. In SON, radiative warming due to decreased α and increased S<sub>dn</sub> offset ~35% of the concurrent evaporative cooling, such that the surface cooling became insignificant.

We also evaluated how the widespread NH greening triggered inter-seasonal T<sub>a</sub> responses. During those seasons when LAI was kept the same in the experiment pair (Table S1), any detected climate anomalies were linked to perturbations of soil and atmospheric conditions arising from LAI variations in the preceding season. At the NH scale, direct surface biophysical processes had poor explanatory power with respect to inter-seasonal changes in T<sub>a</sub> (Fig. S10), signaling weak legacy effects through coupling between soil and surface boundary layer. However, regionally, soil moisture was a key factor connecting MAM LAI changes to JJA T<sub>a</sub> responses, often showing regional warming in areas with soil drying (Figs. S11a, S12a). This occurred because MAM greening enhanced ET, causing a soil moisture deficit that was carried over to JJA; drier soils further enhanced JJA warming by dissipating absorbed solar radiation to a greater extent as sensible rather than latent heat<sup>23</sup>. In other cases, we found similar spatial patterns between trends of 500 hPa geopotential height (z500, an indicator of synoptic circulation changes impacting regional climate) and inter-seasonal ΔT<sub>a</sub>, with negative (positive) z500 trends corresponding to cooling (warming) (Figs. S11, S13). These patterns showed little similarity to greening patterns, suggesting that atmospheric circulation response connected the direct forcing of LAI changes on the local atmosphere to other regions and seasons. Via atmospheric teleconnections, SON NH greening triggered widespread MAM cooling centered in western Asia, as well as local MAM warming centered in northern Canada (Fig. S11g). This MAM cooling matched well with the SON LAI-induced MAM anomalous low, with the western Asia cooling center controlled by a cyclical anomaly structure (Figs. S11g, S13g). By analyzing every repeated experiment with different initial conditions (Methods), we found similar teleconnection patterns in 19 out of 30 ensemble members (Fig. S14), suggesting that this pattern is more likely connected to vegetation greening than internal variability in the climate system. During DJF, an anomalous high prevailed over land forced by greening in all growing seasons (Fig. S13c, f, i), resulting in widespread NH warming (Fig. S11c, f, i), which implies that greening distributed more available solar energy toward the coolest periods of the seasonal cycle. By summarizing the ΔT<sub>a</sub> responses in the space of soil moisture and z500 changes, we found that the inter-seasonal ΔT<sub>a</sub> increases mostly followed a positive gradient of z500 trends (Fig. 5). Therefore, indirect atmospheric processes, such as warm/cold air mass advection and large-scale atmospheric circulation<sup>25,40</sup>, have played a key role in propagating the climatic effects of greening to other seasons.

In summary, this model-based study shows that LAI-T<sub>a</sub> feedbacks through biophysical processes have apparent seasonal inconsistency regarding the signs and magnitudes, resulting from a trade-off between radiative warming and evaporative cooling. This is in contrast to the well-documented carbon-cycle feedbacks, for which NH greening exerts a seasonally consistent negative forcing on T<sub>a</sub> through enhancing plant photosynthetic rates. Our results indicate that, in summer, the overwhelmingly strong ET enhancement forces net surface cooling, whereas in

spring and autumn, radiative warming due to decreased albedo and higher atmospheric water vapor content nullify the evaporative cooling, producing weak and insignificant T<sub>a</sub> changes. Compared with the straightforward response to greening, spring and autumn T<sub>a</sub> changes are impacted to a greater extent by atmospheric circulation anomalies generated by greening in the preceding seasons. We however note potential uncertainties in our model-based LAI-T<sub>a</sub> feedback results, particularly in terms of the inter-seasonal signal, for which solid observational constraints on the atmospheric circulation processes are presently lacking. We encourage other climate model centers to perform similar experiments to refine the exact values of vegetation biophysical feedbacks on seasonal T<sub>a</sub> changes.

Our study achieves a holistic understanding of the seasonal responses of air temperature to NH greening over the growing season. We illustrate complex trade-offs and linkages of vegetation-climate feedbacks between different growing stages, which tend to be concealed within annually aggregated values<sup>13</sup>. This highlights the need to better understand these biophysical processes operating both within and across seasons so that their potential time-lagged climate benefits and/or counterproductive consequences will not be overlooked. The regulatory role of NH greening on seasonal climate also has implications for adaptation planning and decision-making, as greening is now increasingly shaped by human land-use practices such as afforestation and reforestation<sup>41</sup>.

## Materials and methods

**Transient experiments using coupled land-atmosphere models.** We characterized vegetation-climate feedbacks using IPSL-CM, a coupled land-atmosphere GCM developed and maintained by the IPSL modeling community<sup>42</sup>. IPSL-CM has continuously contributed to the fifth and sixth phases of the Coupled Model Intercomparison Project (CMIP), which has informed assessment reports of the Intergovernmental Panel on Climate Change (IPCC)<sup>43,44</sup>. This GCM is composed of the Laboratoire de Météorologie Dynamique general circulation model with zooming capability (LMDZ)<sup>45</sup> and the Organizing Carbon and Hydrology In Dynamic Ecosystems (ORCHIDEE)<sup>46</sup> land surface model. Within the IPSL-CM coupled model, local exchanges of water and energy, as well as large-scale atmospheric circulation changes, are explicitly parameterized. The values of LAI and plant functional types (PFTs) modeled in ORCHIDEE were manually replaced with available satellite observations. Other vegetation-related biophysical parameters such as height, rooting depth, and leaf albedo were determined by prescribed PFT and a set of pre-defined PFT-specific values<sup>46</sup>. The module for calculating the surface energy budget was forced to use prescribed satellite observations, instead of model-simulated values, so observed vegetation changes could alter land-atmospheric energy fluxes and thus modulate T<sub>a</sub>.

Four GCM experiments were conducted to quantify the transient response of T<sub>a</sub> to vegetation greening observed over each of the three growing seasons (MAM, JJA, and SON) during 1982–2014 (Table S1). The four experiments differed in the prescribed seasonal LAI conditions; the SCE run was forced by satellite-observed annually varying monthly LAI maps during 1982–2014, and the three control runs (LAI<sup>MAM</sup><sub>CTL</sub>, LAI<sup>IIA</sup><sub>CTL</sub>, and LAI<sup>SON</sup><sub>CTL</sub>, corresponding to the three growing seasons) used the same LAI maps except for the season of interest, for which monthly varying climatological LAI maps for 1982–2014 were used. In each growing season, the difference between the experiment pair, with yearly varying and fixed present-season LAI (SCE-LAI<sup>MAM</sup><sub>CTL</sub> for MAM, SCE-LAI<sup>IIA</sup><sub>CTL</sub> for JJA, and SCE-LAI<sup>SON</sup><sub>CTL</sub> for SON), distinguished the biophysical feedbacks of observed seasonal LAI changes on the surface energy budget and T<sub>a</sub> of both the present and subsequent seasons.

The satellite-observed LAI maps used for the experiments were obtained from the Global Inventory Monitoring and Modeling Studies (GIMMS) LAI3g, which are available at an 8-km spatial resolution and 15-day intervals for 1982–2014 (ref. 36). This global LAI product was generated from GIMMS AVHRR Normalized Difference Vegetation Index 3 g data using an artificial neural network algorithm<sup>36</sup>. The GIMMS LAI3g data replace possible snow- or cloud-contaminated LAI values with an average seasonal profile or interpolated values<sup>36</sup>. The high quality of this data product has been demonstrated through direct comparisons with field measurements<sup>36</sup>. This LAI product has been widely used to study the seasonal evolution of vegetation activity and its interactions with climate<sup>23,37,47,48</sup>. Moreover, other AVHRR-based LAI products (e.g., GLOBMAP<sup>49</sup> and GLASS<sup>50</sup>) are available in the literature but have a considerable discrepancy regarding long-term LAI trends. Among these products, GIMMS LAI shows intermediate rates of greening seasonally<sup>5</sup>, thus offering a moderate estimate for LAI-climate feedback.
As satellite-based land cover maps are not available for the entire study period, we used the static 0.5-degree Olson land cover map generated by the ORCHIDEE
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Biophysical processes underlying inter-seasonal temperature responses to greening. Distribution of the ΔT<sub>a</sub> (air temperature) trend in the space of the ΔSM (soil moisture) trend and the Δz500 (500 hPa geopotential height) trend, all induced by MAM (**a**, **d**, **g**), JJA (**b**, **e**, **h**), and SON (**c**, **f**, **i**) leaf area index (LAI) changes. Marginal curves indicate the partial dependence of the ΔT<sub>a</sub> trend on the ΔSM trend (top) and the Δz500 trend (left). MAM March-April-May, JJA June-July-August, SON September-October-November, DJF December-January-February.</figcaption></figure>

community (https://forge.ipsl.jussieu.fr/orchidee/wiki/Documentation/Ancillary). The Olson map originally contains 94 land categories at 5 km resolution⁵¹, which were converted to the 12 PFTs and bare soil using a look-up table approach. The use of static PFT maps may to some degree bias the model results for regions where the LAI trends are dominated by human activities (e.g., agricultural intensification, plantation, deforestation), but will not have a hemispheric effect since the seasonal greening is contributed mainly by rising CO₂ and climate change⁵. To meet the model input requirements, monthly LAI maps were first calculated as maximum values of 15-day LAI composites and then aggregated to 0.5° × 0.5° grid cells. Next, LAI maps at the PFT level were derived by weighting the satellite-observed LAI values by the fractional vegetation cover of each PFT within each grid cell.

Apart from seasonal LAI differences, all experiments were forced identically with realistic SST, SIC, and atmospheric CO₂ concentrations observed from 1982–2014. The inclusion of observed surface boundary conditions allowed the GCM to better reproduce observed interannual and longer-term variations in global Tₐ. Therefore, the separated vegetation-temperature feedbacks encompassed the indirect effects of vegetation–atmosphere–ocean interactions and CO₂ radiative and physiological forcings, which are, in reality, embedded in observed temperature changes. Monthly SST and SIC maps were obtained from the Atmospheric Model Intercomparison Project (AMIP; http://www-pcmdi.llnl.gov/projects/amip/) with 1° × 1° spatial resolution. Observed global annual atmospheric CO₂ concentrations were derived from those used for transient simulations in the “Trends and drivers of the regional scale sources and sink of carbon dioxide” (TRENDY) project (http://dgym.ceh.ac.uk/node/9).

Figure S1 shows the steps for performing model simulations. During model spin-up, we performed 30-year undisturbed runs forced by climatological LAI, SST, SIC, and atmospheric CO₂ to ensure that soil moisture and temperature reached an equilibrium with the climate conditions. Next, we extended the undisturbed runs for

another 30 years to generate a 30-member initial-condition ensemble, with outputs of each year for providing an alternative initial condition. This approach is effective for minimizing uncertainties in modeled climate responses to perturbations related to internal climate variability¹³,⁵². Starting with these different initial conditions, we performed the four transient simulations using different LAI maps (Table S1) from January 1982 to December 2014 as input. We calculated the LAI-induced Tₐ changes (i.e., ΔTₐ) as the difference between SCE and the simulation with fixed seasonal LAI (LAI<sup>MAM</sup><sub>CTL</sub>, LAI<sup>JJA</sup><sub>CTL</sub>, and LAI<sup>SON</sup><sub>CTL</sub>), both using the 30-member ensemble mean of the transient simulations. For each seasonal experiment, intra-seasonal Tₐ changes were quantified in the same season of fixed LAI, and inter-seasonal Tₐ changes (carry-over effects) were quantified for the remaining seasons. Only vegetated northern areas with climatological LAI ≥ 0.1 were considered in this analysis.

We also conducted numerical experiments using another coupled land–atmosphere model, CAM-CLM, to verify the robustness of our results. We used the model that couples the Community Atmosphere Model, version 6 (CAM6) and the Community Land Model, version 5 (CLM5)⁵³, with a spatial resolution of 1.9° × 2.5°. Four 100-year-long equilibrium experiments were conducted with CAM-CLM using different LAI inputs: the LAI<sub>2010s</sub> was prescribed with seasonally varying LAI averaged for the 2010–2014 period; and the three seasonal simulations (LAI<sup>MAM</sup><sub>1980s</sub>, LAI<sup>JJA</sup><sub>1980s</sub>, and LAI<sup>SON</sup><sub>1980s</sub>) used the same LAI maps except for the focused season for which the average LAI for 1982–1986 was used (Table S2). The LAI-temperature feedbacks were quantified as the difference between the LAI<sub>2010s</sub> and the seasonal simulations, averaged over the last 30 years of the 100-year-long runs.

**Observation-based T<sub>a</sub> used for model evaluation.** We used observation-based T<sub>a</sub> data to evaluate the performance of the IPSL-CM GCM in simulating NH warming patterns and their seasonality. Air temperature data were derived from the
Princeton Global Forcing (PGF) product developed by the Terrestrial Hydrology Research Group at Princeton University<sup>54,55</sup>. This forcing dataset provides near-surface meteorology data generated from the National Centers for Environmental Prediction-National Center for Atmospheric Research (NCEP-NCAR) reanalysis, which was bias-corrected based on the Climatic Research Unit (CRU) temperature<sup>54</sup>. The PGF temperature data are available at 3-h and 1° × 1° resolution covering the period 1948–2014.

**Analytical decomposition of LAI-induced T<sub>a</sub> changes.** We employed the analytical approaches of ref. ¹³ and ref. ³⁷ to disentangle the contributions of key surface biophysical parameters to the modeled LAI-induced seasonal changes in surface radiative forcing. The considered biophysical processes and parameters included λE (latent heat flux, which cools the land surface by transferring absorbed energy to the troposphere), α (surface albedo, the ratio of reflected shortwave radiation to incoming shortwave radiation), S<sub>dn</sub> (surface downward shortwave radiation, controlled by cloud and atmospheric water vapor content), ε<sub>a</sub> (air emissivity, which reflects the effectiveness of air in emitting and absorbing long-wave radiation from the surface), and r<sub>a</sub> (aerodynamic resistance, which affects the efficiency of turbulent transfer of heat and water).

The decomposition was based inherently on the surface energy balance. Briefly, vegetation changes affect the surface energy budget and alter the T<sub>s</sub>, further warming or cooling the local near-surface air. The surface energy balance equation is written as:

$$S_{\text{net}} + L_{\text{net}} = H + \lambda E + G \quad (1)$$

where S<sub>net</sub>, L<sub>net</sub>, H, and G are the surface net shortwave radiation, surface net longwave radiation, sensible heat flux, and ground heat flux, respectively. The left-hand side of Eq. (1) represents the net solar and thermal radiation absorbed by the surface (i.e., radiative processes), and the right-hand side represents the energy dissipation at the surface (i.e., non-radiative processes). G can be neglected on seasonal and interannual timescales, and the other terms can be expressed as follows:

$$S_{\text{net}} = S_{\text{dn}} - S_{\text{up}} = S_{\text{dn}}(1 - \alpha) \quad (2)$$

$$L_{\text{net}} = L_{\text{dn}} - L_{\text{up}} = \epsilon_s \sigma (\epsilon_a T_a^4 - T_s^4) \quad (3)$$

$$H = \rho C_p \frac{T_s - T_a}{r_a} \quad (4)$$

where S<sub>dn</sub> and S<sub>up</sub> are the surface downward and upward shortwave radiation, respectively, and L<sub>dn</sub> and L<sub>up</sub> are the surface downward and upward longwave radiation, respectively. Parameters ε<sub>s</sub>, σ, ρ, and C<sub>p</sub> are the land surface emissivity (0.95), Stephan–Boltzmann constant (5.67 × 10<sup>-8</sup> W m<sup>-2</sup> K<sup>-4</sup>), air density (1.205 kg m<sup>-3</sup>), and specific heat capacity of the air at constant pressure (1013 J kg<sup>-1</sup> K<sup>-1</sup>), respectively.

By combining Eqs. (2–4) and differentiating Eq. (1) with respect to T<sub>s</sub>, ΔT<sub>s</sub> can be derived as follows:

$$\Delta T_s = \frac{1}{f_s} \left( -S_{\text{dn}}\Delta\alpha + (1-\alpha)\Delta S_{\text{dn}} - \Delta\lambda E + \frac{\rho C_p}{r_a^2} (T_s - T_a)\Delta r_a + \epsilon_s \sigma T_a^4 \Delta\epsilon_a \right) + \left( \frac{\rho C_p}{r_a} + 4\epsilon_s \sigma \epsilon_a T_a^3 \right) \Delta T_a / \left( \frac{\rho C_p}{r_a} + 4\epsilon_s \sigma T_a^3 \right) \quad (5)$$

$$f_s = \frac{\rho C_p}{r_a} + 4\epsilon_s \sigma T_s^3 \quad (6)$$

where f<sub>s</sub> is the surface energy redistribution factor, which represents changes in T<sub>a</sub> in response to radiative forcing of 1 W m<sup>-2</sup> at the land surface. The first term on the right-hand side of Eq. (5) denotes LAI-induced T<sub>s</sub> changes due to altered surface radiative forcing (ΔT<sup>rad</sup><sub>s</sub> = f<sub>s</sub><sup>-1</sup>(-S<sub>dn</sub>Δα + (1-α)ΔS<sub>dn</sub> - ΔλE + ρC<sub>p</sub>(T<sub>s</sub> - T<sub>a</sub>)r<sub>a</sub><sup>-2</sup>Δr<sub>a</sub> + ε<sub>s</sub>σT<sub>a</sub><sup>4</sup>Δε<sub>a</sub>), and the second term demonstrates that LAI-induced T<sub>s</sub> changes can be almost equivalently converted to T<sub>a</sub> changes (ΔT<sup>rad</sup><sub>a</sub> = (ρC<sub>p</sub>r<sub>a</sub><sup>-1</sup> + 4ε<sub>s</sub>σT<sub>s</sub><sup>3</sup>)ΔT<sup>rad</sup><sub>s</sub> / (ρC<sub>p</sub>r<sub>a</sub><sup>-1</sup> + 4ε<sub>s</sub>σT<sub>a</sub><sup>3</sup>)). Local ΔT<sub>a</sub> changes unexplained by direct surface-atmospheric energy exchange (ΔT<sup>rad</sup><sub>a</sub>) can be attributed to air mass advection or changes in atmospheric circulation (ΔT<sup>cir</sup><sub>a</sub>), as follows:

$$\Delta T_a = \frac{1}{f_a} \left( -S_{\text{dn}}\Delta\alpha + (1 - \alpha)\Delta S_{\text{dn}} - \Delta\lambda E + \frac{\rho C_p}{r_a^2} (T_s - T_a)\Delta r_a + \epsilon_s \sigma T_a^4 \Delta\epsilon_a \right) + \Delta T_a^{\text{cir}} \quad (7)$$

where f<sub>a</sub> is the surface energy redistribution factor for T<sub>a</sub>, given by

$$f_a = \frac{\rho C_p}{r_a} + 4\epsilon_s \sigma \epsilon_a T_s^3 \quad (8)$$

Similar to previous assessments<sup>13,37</sup>, climatic effects of vegetation changes associated with different biophysical processes were assessed with equivalent surface radiative forcing, to ensure direct comparability among them. In Eq. (7), -S<sub>dn</sub>Δα, (1-α)ΔS<sub>dn</sub>, -ΔλE, ρC<sub>p</sub>r<sub>a</sub><sup>-2</sup>(T<sub>s</sub> - T<sub>a</sub>)Δr<sub>a</sub>, and ε<sub>s</sub>σT<sub>a</sub><sup>4</sup>Δε<sub>a</sub> represent surface

radiative forcing caused by vegetation-induced changes of α, S<sub>dn</sub>, λE, r<sub>a</sub>, and ε<sub>a</sub>, respectively. The sum of -S<sub>dn</sub>Δα, (1-α)ΔS<sub>dn</sub>, -ΔλE and ε<sub>s</sub>σT<sub>a</sub><sup>4</sup>Δε<sub>a</sub> represents radiative forcing associated with radiative processes, and the sum of -ΔλE and ρC<sub>p</sub>r<sub>a</sub><sup>-2</sup>(T<sub>s</sub> - T<sub>a</sub>)Δr<sub>a</sub> represents radiative forcing associated with non-radiative processes.

## Data availability

All observation and model data that support the findings of this study are available as follows. The AVHRR GIMMS LAI3g data were available at https://cliveg.bu.edu/modismirsl/laI3g-fpar3g.html. The climatic variables from the Princeton Global Meteorological Forcing (PGF) v3.0 data were available at https://hydrology.princeton.edu/data.pgf.php/.

## Code availability

Code and documentation for the IPSL model are publicly available at https://cmc.ipsl.fr/ipsl-climate-models/ipsl-cm5/. Description and code relevant to CAM6 and CLM5 are freely available online at https://www.cesm.ucar.edu/models/cesm2. Model outputs were processed using the software Matlab 2019a. The MATLAB scripts are available from the corresponding author upon request.

Received: 30 October 2021; Accepted: 20 June 2022; Published online: 07 July 2022

## References

1. Keenan, T. F. et al. Net carbon uptake has increased through warming-induced changes in temperate forest phenology. *Nat. Clim. Change* **4**, 598–604 (2014).

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

3. Liu, Q. et al. Delayed autumn phenology in the Northern Hemisphere is related to change in both climate and spring phenology. *Glob. Change Biol.* **22**, 3702–3711 (2016).

4. Richardson, A. D. et al. Climate change, phenology, and phenological control of vegetation feedbacks to the climate system. *Agric. Forest Meteorol.* **169**, 156–173 (2013).

5. Zhu, Z. et al. Attribution of seasonal leaf area index trends in the northern latitudes with “optimally” integrated ecosystem models. *Glob. Change Biol.* **23**, 4798–4813 (2017).

6. Piao, S. et al. Characteristics, drivers and feedbacks of global greening. *Nat. Rev. Earth Environ.* **1**, 14–27 (2020).

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

8. Lian, X. et al. Seasonal biological carryover dominates northern vegetation growth. *Nat. Commun.* **12**, 983 (2021).

9. Peñuelas, J., Rutishauser, T. & Filella, I. Ecology phenology feedbacks on climate change. *Science* **324**, 887–888 (2009).

10. Barichivich, J., Briffa, K. R., Osborn, T. J., Melvin, T. M. & Caesar, J. Thermal growing season and timing of biospheric carbon uptake across the Northern Hemisphere. *Glob. Biogeochem. Cycles* **26**, GB4015 (2012).

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

12. Schwartz, M. D. & Karl, T. R. Spring phenology: nature's experiment to detect the effect of “green-up” on surface maximum temperatures. *Mon. Weather Rev.* **118**, 883–890 (1990).

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

14. Lee, X. et al. Observed increase in local cooling effect of deforestation at higher latitudes. *Nature* **479**, 384–387 (2011).

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

16. Hogg, E. H., Price, D. & Black, T. Postulated feedbacks of deciduous forest phenology on seasonal climate patterns in the western Canadian interior. *J. Clim.* **13**, 4229–4243 (2000).

17. Moore, K. E. et al. Seasonal variation in radiative and turbulent exchange at a deciduous forest in central Massachusetts. *J. Appl. Meteorol. Rev.* **35**, 122–134 (1996).

18. Betts, R. A. Offset of the potential carbon sink from boreal forestation by decreases in surface albedo. *Nature* **408**, 187–190 (2000).

19. Abe, M., Takata, K., Kawamiya, M. & Watanabe, S. Vegetation masking effect on future warming and snow albedo feedback in a boreal forest region of
northern Eurasia according to MIROC-ESM. *J. Geophys. Res. Atmos.* **122**, 9245–9261 (2017).

20. Loranty, M. M., Berner, L. T., Goetz, S. J., Jin, Y. & Randerson, J. T. Vegetation controls on northern high latitude snow-albedo feedback: observations and CMIP5 model simulations. *Glob. Change Biol.* **20**, 594–606 (2014).

21. Xu, X., Riley, W. J., Koven, C. D., Jia, G. & Zhang, X. Earlier leaf-out warms air in the north. *Nat. Clim. Change* **10**, 370–375 (2020).

22. Swann, A. L., Fung, I. Y., Levis, S., Bonan, G. B. & Doney, S. C. Changes in Arctic vegetation amplify high-latitude warming through the greenhouse effect. *Proc. Natl Acad. Sci. USA* **107**, 1295–1300 (2010).

23. Lian, X. et al. Summer soil drying exacerbated by earlier spring greening of northern vegetation. *Sci. Adv.* **6**, eaax0255 (2020).

24. Wolf, S. et al. Warm spring reduced carbon cycle impact of the 2012 US summer drought. *Proc. Natl Acad. Sci. USA* **113**, 5880–5885 (2016).

25. Winckler, J., Lejeune, Q., Reick, C. H. & Pongratz, J. Nonlocal effects dominate the global mean surface temperature response to the biogeophysical effects of deforestation. *Geophys. Res. Lett.* **46**, 745–755 (2019).

26. Xia, J. et al. Terrestrial carbon cycle affected by non-uniform climate warming. *Nat. Geosci.* **7**, 173–180 (2014).

27. Stine, A. R., Huybers, P. & Fung, I. Y. Changes in the phase of the annual cycle of surface temperature. *Nature* **457**, 435–440 (2009).

28. Mann, M. E. & Park, J. Greenhouse warming and changes in the seasonal cycle of temperature: Model versus observations. *Geophys. Res. Lett.* **23**, 1111–1114 (1996).

29. Alkama, R. & Cescatti, A. Biophysical climate impacts of recent changes in global forest cover. *Science* **351**, 600–604 (2016).

30. Peng, S. et al. Afforestation in China cools local land surface temperature. *Proc. Natl Acad. Sci. USA* **111**, 2915–2919 (2014).

31. Chen, L. & Dirmeyer, P. A. Reconciling the disagreement between observed and simulated temperature responses to deforestation. *Nat. Commun.* **11**, 202 (2020).

32. Devaraju, N., de Noblet-Ducoudré, N., Quesada, B. & Bala, G. Quantifying the relative importance of direct and indirect biophysical effects of deforestation on surface temperature and teleconnections. *J. Clim.* **31**, 3811–3829 (2018).

33. Findell, K. L. et al. The impact of anthropogenic land use and land cover change on regional climate extremes. *Nat. Commun.* **8**, 989 (2017).

34. Park, H. & Jeong, S. Leaf area index in Earth system models: how the key variable of vegetation seasonality works in climate projections. *Environ. Res. Lett.* **16**, 034027 (2021).

35. Marti, O. et al. Key features of the IPSL ocean atmosphere model and its sensitivity to atmospheric resolution. *Clim. Dyn.* **34**, 1–26 (2010).

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

37. Li, Y., Piao, S., Chen, A., Ciais, P. & Li, L. Z. X. Local and tele-connected temperature effects of afforestation and vegetation greening in China. *Natl Sci. Rev.* **7**, 897–912 (2019).

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

39. Bateni, S. M. & Entekhabi, D. Relative efficiency of land surface energy balance components. *Water Resour. Res.* **48**, 4510 (2012).

40. Swann, A. L. S. et al. Continental-scale consequences of tree die-offs in North America: identifying where forest loss matters most. *Environ. Res. Lett.* **13**, 055014 (2018).

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

42. Dufresne, J. L. et al. Climate change projections using the IPSL-CM5 Earth System Model: from CMIP3 to CMIP5. *Clim. Dyn.* **40**, 2123–2165 (2013).

43. Taylor, K. E., Stouffer, R. J. & Meehl, G. A. An overview of CMIP5 and the experiment design. *Bull. Am. Meteorol. Soc.* **93**, 485–498 (2012).

44. Eyring, V. et al. Overview of the coupled model intercomparison project phase 6 (CMIP6) experimental design and organization. *Geosci. Model Dev.* **9**, 1937–1958 (2016).

45. Hourdin, F. et al. The LMDZ4 general circulation model: climate performance and sensitivity to parametrized physics with emphasis on tropical convection. *Clim. Dyn.* **27**, 787–813 (2006).

46. Krinner, G. et al. A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system. *Glob. Biogeochem. Cycles* **19**, GB1015 (2005).

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

48. Winkler, A. J., Myneni, R. B., Alexandrov, G. A. & Brovkin, V. Earth system models underestimate carbon fixation by plants in the high latitudes. *Nat. Commun.* **10**, 885 (2019).

49. Liu, R. & Liu, Y. Generation of new cloud masks from MODIS land surface reflectance products. *Remote Sens. Environ.* **133**, 21–37 (2013).

50. Liang, S. et al. The global land surface satellite (GLASS) product suite. *Bull. Am. Meteorol. Soc.* **102**, E323–E337 (2021).

51. Olson, D. M. et al. Terrestrial ecoregions of the world: a new map of life on Earth. *BioScience* **51**, 933–938 (2001).

52. Deser, C. et al. Insights from Earth system model initial-condition large ensembles and future prospects. *Nat. Clim. Change* **10**, 277–286 (2020).

53. Danabasoglu, G. et al. The Community Earth System Model version 2 (CESM2). *J. Adv. Model. Earth Syst.* **12**, e2019MS001916 (2020).

54. Sheffield, J., Goteti, G. & Wood, E. F. Development of a 50-year high-resolution global dataset of meteorological forcings for land surface modeling. *J. Clim.* **19**, 3088–3111 (2006).

55. Bright, R. M. et al. Local temperature response to land cover and management change driven by non-radiative processes. *Nat. Clim. Change* **7**, 296–302 (2017).

## Acknowledgements

This study was supported by the Second Tibetan Plateau Scientific Expedition and Research (STEP) program (grant no. 2019QZKK0208), the National Natural Science Foundation of China (41988101), and the Xplorer Prize. S.J. was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT) (NRF-2019R1A2C3002868). J.P. was supported by the Spanish Government grant PID2019-110521GB-100, the Fundación Ramón Areces grant ELEMENTAL-CLIMATE, and the Catalan Government grants SGR2017-1005 and AGAUR-2020PANDE00117.

## Author contributions

X.L., S.P., and S.J. designed the research; X.L. performed analyses, drafted the figures, and wrote the first draft of the manuscript; C.-E.P., H.X., L.Z.X.L., T.W., P.G., and J.P. contributed to the interpretation of the results and to the text.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-022-31671-z.

**Correspondence** and requests for materials should be addressed to Sujong Jeong or Shilong Piao.

**Peer review information** *Nature Communications* thanks Heping Liu, and the other, anonymous, reviewers for their contribution to the peer review of this work.

**Reprints and permission information** is available at http://www.nature.com/reprints

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

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022