Main

Increasing frequency and severity of extreme droughts can critically damage the stability and health of terrestrial ecosystems and weaken the land carbon sink1,2,3. The time required for an ecosystem to recover from severe drought and return to its predrought functional state (hereafter termed drought recovery) is an important component of ecosystem resilience4,5,6,7. Drought recovery from a single extreme drought event, within the growing season, has been widely studied7,8,9,10. However, interactions between vegetation phenology (the start and end of the growing season) and drought timing and intensity have been understudied and remain largely unknown (Fig. 1)11,12,13. Vegetation phenology can have crucial effects on drought recovery through biophysical mechanisms. For instance, changes in spring phenology drive evapotranspiration and thus alter the amount of water available to support vegetation growth in subsequent seasons14. Vegetation phenology can also impact drought recovery through biological processes, whereby changes in spring phenology directly influence the health and physiology of vegetation, which mediates the ability of vegetation to recover from drought15.

Fig. 1: Schematic diagram to quantify drought recovery and the governing processes.
figure 1

Two different trajectories of drought recovery (time required for an ecosystem to recover to its predrought functional state from extreme drought events) were specified depending on the interacting effects of vegetation phenology, drought timing and other bioclimatic and soil factors. a, Drought recovery can occur within a single growing season (RSGS) or extends to the subsequent growing season (drought recovery within multiple growing seasons) (RMGS). Under RMGS, the drought recovery is suspended by vegetation dormancy and is strongly dependent on climate conditions during the dormancy period. Both the climate and vegetation growth conditions in the 6 months preceding extreme drought events (predrought) were considered. Drought response lag was defined as the duration (months) between the onset of extreme drought event and the deepest suppression of vegetation growth induced from the drought events. b, Vegetation phenology effects on drought recovery through both biophysical and biological effects, interacted with multiple bioclimatic and soil factors as well as the drought sensitivity. SOSdrought_year and SOSsubsequent_year are the start of growing season in the drought year and the subsequent year, respectively, and EOSdrought_year is the end of growing season in the drought year. The progression of seasons is from left to right in a and b. Only single extreme drought events are considered in this study (Methods).

Vegetation phenology effects on drought recovery are probably pervasive in the mid- and high-latitudinal Northern Hemisphere, where vegetation growth is characterized by a marked phenological cycle12,16,17. Drought recovery is sensitive to and dependent on multiple complex and interacting factors including vegetation phenology, drought timing and anomalies in seasonal bioclimatic factors. These factors impact the overall timing of drought recovery which can occur relatively quickly over a single growing season (RSGS) or persist to subsequent growing seasons (within multiple growing seasons) (RMGS) (Methods). Vegetation phenology interacts with drought timing and thus determines subsequent vegetation growth11,14,15 (Fig. 1). On the one hand, spring phenology strongly governs vegetation growth and soil water consumption in the growing season11,12,13,18, thus influencing drought recovery. On the other hand, changes in temperature and snow accumulation (snow water equivalent (SWE)), during the non-growing season can directly determine subsequent spring phenology through biological processes19,20, thus influencing drought recovery (Fig. 1). Vegetation phenology is already being rapidly impacted by climate warming12,21 and factors such as shifting drought seasonality are already reshaping the vulnerability of ecosystems to intensified extreme drought events to a currently unknown degree. Quantifying these complex feedbacks requires a systematic understanding of the interactions between vegetation phenology, drought timing, bioclimatic factors and drought recovery (Fig. 1).

Here, we quantify the vegetation phenology effects on drought recovery over the mid- and high-latitudinal Northern Hemisphere—a region with pronounced seasonality and susceptibility to drought21. Changes in vegetation greenness caused by droughts are identified by satellite observations of normalized difference vegetation index (NDVI; 1982–2015)22, complemented by two additional independent ecophysiological proxies of vegetation dynamics—contiguous sun-induced fluorescence (CSIF; 2001–2015)23,24 and microwave-based vegetation optical depth (VOD; 2003–2015)25 (Methods). Extreme drought events are defined using the multiscalar standardized precipitation evapotranspiration index (SPEI26; Methods). The main objectives are: (1) to identify the divergent trajectories in drought recovery (Fig. 1a) and (2) to quantify the impacts of vegetation phenology on drought recovery (Fig. 1b). We further quantify the role of multiple bioclimatic factors, including pre- and postdrought climatic conditions, preceding growth status, ecophysiological traits and soil properties in mediating the vegetation phenology effects on drought recovery. We test the overarching hypothesis that spring phenology has critical effects on drought recovery that depends on the timing of drought. Specifically, we hypothesize that: (1) an earlier spring phenology in the drought year will drive longer drought recovery across vegetation types due to the dominant role of seasonal biophysical feedbacks and (2) a delayed spring phenology in the subsequent year will postpone vegetation growth, resulting in longer drought recovery across vegetation types due to the dominant role of biological processes (Fig. 1).

Drought recovery pattern

We first quantified and compared the spatial pattern in the two trajectories of drought recovery: RSGS and RMGS (Fig. 2 and Supplementary Fig. 1; Methods). Unexpectedly, our analysis revealed that in ~50% (ranging from 45% to 51% among evergreen forests, deciduous forests, shrubs and grasses) of early-growing season extreme drought events, NDVI-based vegetation greenness did not fully recover to the predrought condition within a single growing season (Fig. 2f). This percentage increased to more than 60% and 80% when extreme drought events occurred in the mid- and late-growing season, respectively (Fig. 2f). Longer drought recovery was found in central North America, the Mediterranean and central Eurasia under both RSGS and RMGS trajectories (Fig. 2a–c). Across these regions, mean growing season NDVI, CSIF and VOD exhibited consistent and significantly positive correlation with mean growing season SPEI over a variety of timescales (Extended Data Fig. 1). By contrast, in the northern latitudes (>50° N), in 84% of extreme drought events, NDVI-based vegetation greenness recovered in <3 months (Fig. 2a). Our analysis further revealed that mean drought recovery time decreased with latitude and aridity (Extended Data Fig. 2 and Supplementary Fig. 2). Consistent spatial patterns in drought recovery were also obtained from different definitions of extreme drought events using SPEI aggregated over timescales of 1, 6, 9 or 12 months (Supplementary Fig. 3). Additionally, consistent patterns in drought recovery were observed from the independent vegetation growth proxies of CSIF and VOD (Extended Data Fig. 3).

Fig. 2: Spatial patterns in drought recovery in mid- and high-latitudinal Northern Hemisphere.
figure 2

ac, The spatial pattern of mean drought recovery (months) averaging both RSGS and RMGS (a), drought recovery under RSGS (b) and drought recovery under RMGS (c). d,e, Drought recovery in different combinations of drought timing with different vegetation types and all vegetated area under RSGS (d) and RMGS (e) scenarios. f, Percentage (%) of extreme drought events that vegetation growth fails to fully recover within a single growing season for different combinations of drought timing and vegetation types. The drought recovery from late-growing season drought under RSGS was not considered to prevent the inaccuracy possibly caused by insufficient training dataset samples (d). Error bars in d and e indicate the 95% CIs of the mean drought recovery. Different capital letters in d and e indicate statistical significance (P < 0.05) in drought recovery among vegetation types with the same drought timing, while different lower-case letters indicate the same vegetation type but different drought timing by two-sided Student’s t-test. We grouped extreme drought events into early-growing season drought (early-GS), mid-growing season drought (mid-GS) and late-growing season drought (late-GS).

Source data

Overall, drought recovery following early-growing season droughts was on average longer than that following mid- and late-growing season droughts. Under RSGS, drought recovery following early-growing season droughts was 0.4 months longer than that of mid-growing season droughts (P < 0.05) (Fig. 2d). Under RMGS, drought recovery following early-growing season droughts was 0.6 and 1.8 months longer than that of mid- and late-growing season droughts, respectively (P < 0.05) (Fig. 2e). Across vegetation types, drought recovery was the longest in deciduous forests compared to evergreen forests, grasses and shrubs in early-, mid- and late-growing season droughts under both RSGS and RMGS (Fig. 2d,e).

Drivers of drought recovery

We tested the hypothesis that spring phenology had significant effects on drought recovery by quantifying the drivers governing drought recovery among different vegetation types under both RSGS and RMGS. The random forest (RF) models were applied to attribute drought recovery to NDVI-derived vegetation canopy phenology, drought sensitivity (the dependency of mean growing season NDVI on mean growing season drought condition during 1982–2015) and multiple bioclimatic and soil factors (Methods). Separate RF models were built for each of the eight groups including two recovery modes (RSGS and RMGS) and four vegetation types (evergreen forests, deciduous forests, shrubs and grasses). The RF models performed well in capturing drought recovery, with better performance for RMGS (R2 ranged from 0.86 to 0.90 across vegetation types) than for RSGS (R2 ranged from 0.71 to 0.75).

The analyses identified significant but differentiated phenological effects on drought recovery among different vegetation types and between RSGS and RMGS. For RSGS, in supporting the first hypothesis, spring phenology was found to strongly impact drought recovery (Fig. 3e,g). However, the effect of spring phenology on drought recovery displayed a bimodal response: both earlier and delayed spring phenology in drought year postponed drought recovery within a single growing season for all vegetation types (Fig. 3e). We further tested whether this bimodal pattern was caused by the interaction between spring phenology and drought timing. Alternative RF models were constructed for each vegetation type under RSGS, while accounting for different drought timings. The analyses showed that earlier spring phenology during the drought year led to shorter drought recovery when extreme drought occurred in the early-growing season but it led to longer drought recovery when extreme drought occurred during the mid-growing season (Extended Data Fig. 4).

Fig. 3: Response function for the drought recovery within a single growing season.
figure 3

af, Response function along each independent variable (x axis), including total precipitation during postdrought period (Precippost) (a), mean vapour pressure deficit during postdrought period (VPDpost) (b), mean NDVI in the 6 months preceding the extreme drought events (NDVIpre) (c), mean annual temperature (MAT, °C) (d), the start of growing season of drought year (SOSdrought_year, d) (e) and mean VPD during the drought response lag period (VPDlag) (f). The y axis in af shows the drought recovery (months). Shaded areas of different colours in af are the 95% CIs for different vegetation types. The bootstrap method was used to compute the average over all the predictions of drought recovery and the 95% CI. g, Normalized variable importance (unitless) for the drought recovery of different vegetation types. Four kinds of vegetation types were considered and compared, including evergreen forests, deciduous forests, shrubs and grasses (Methods). Precippre, total precipitation in the 6 months preceding the extreme drought events; VPDpre, mean VPD in the 6 months preceding the extreme drought events; Preciplag, total precipitation during the drought response lag period; MAP, mean annual precipitation; PrecipCV, coefficient of variation of mean annual precipitation; and TempCV, coefficient of variation of mean annual temperature. Numbers in the legend indicate the sampling sizes. The z-scores were calculated for Precippost, VPDpost, NDVIpre, VPDlag, Precippre, VPDpre and Preciplag before the RF analysis.

Source data

We next tested the second hypothesis that a delayed spring phenology in the subsequent year will postpone drought recovery. For RMGS, in support of the second hypothesis, delayed spring phenology in the subsequent year consistently delayed drought recovery, while earlier spring phenology shortened drought recovery, with importance of spring phenology ranging from 46% to 58% across vegetation types (Fig. 4a,i). Climate conditions during the dormancy period also played important roles for RMGS but were less important than spring phenology (Fig. 4e,f,i). A negative temperature anomaly during the dormancy period tended to delay drought recovery for all vegetation types (Fig. 4f). Interestingly, both negative and positive anomalies in SWE led to longer drought recovery, with this effect more prominent in deciduous forests and shrubs (with importance of 20% and 19%, respectively) (Fig. 4i). We further found that a longer dormant season was associated with shorter drought recovery (Fig. 4d). Negative anomalies in precipitation (with importance of 52% and 34%, respectively) and positive anomalies in VPD (with importance of 39% and 23%, respectively) during the postdrought period delayed drought recovery under both RSGS and RMGS (Figs. 3a,b and 4b,c). Higher mean annual temperature led to longer drought recovery for both RSGS and RMGS, irrespective of vegetation types (Figs. 3d and 4h). However, total precipitation in the 6 months preceding extreme drought events and the multiyear mean water balance and precipitation only exerted marginal effects on drought recovery (Figs. 3g and 4i). We did not find important effects of either a drought response lag (the time lag between onset of extreme drought and the maximum vegetation growth reduction) or precipitation and VPD during this lag period, on drought recovery for either RSGS or RMGS. The interannual variability of mean annual temperature and mean annual precipitation also showed no significant effects on drought recovery for either RSGS or RMGS (Figs. 3g and 4i).

Fig. 4: Response function for the drought recovery extended to subsequent year.
figure 4

ah, Response function along each independent variable (x axis), including the start of growing season of the subsequent year (SOSsubsequent_year, d) (a), total precipitation during postdrought period (Precippost) (b), mean vapour pressure deficit during postdrought period (VPDpost) (c), dormant length (d) (d), SWE during dormant period (SWEdorm) (e), mean temperature during dormant period (Tempdorm) (f), mean NDVI in 6 months preceding the extreme drought events (NDVIpre) (g) and mean annual temperature in period 1982–2015 (MAT, °C) (h). Shaded areas of different colours in ah are the 95% CIs for different vegetation types. The bootstrap method was used to compute the average over all the predictions of drought recovery and the 95% CI. The y axis in ah shows the drought recovery (months). i, Normalized variable importance (unitless) for the drought recovery of different vegetation types. Four kinds of vegetation types were considered and compared, including evergreen forests, deciduous forests, shrubs and grasses (Methods). Precippre, total precipitation in 6 months preceding the extreme drought events; SOSdrought_year, the start of growing season of drought year; MAP, mean annual precipitation; PrecipCV, coefficient of variation of mean annual precipitation; TempCV, coefficient of variation of mean annual temperature; VPDpre, mean VPD in 6 months preceding the extreme drought events; Preciplag, total precipitation during the drought response lag period; and VPDlag, mean VPD during the drought response lag period. Numbers in the legend indicate the sampling size. z-scores were calculated for Precippost, VPDpost, SWEdorm, Tempdorm, NDVIpre, Precippre, VPDpre, Preciplag and VPDlag before the RF analysis.

Source data

Finally, we found that the NDVIpre (mean NDVI anomaly in the 6 months within growing season preceding extreme drought) exerted important effects on drought recovery, with more prominent effects in RSGS (with importance of 25–31% across vegetation types). For instance, negative anomalies in NDVIpre tended to delay drought recovery under both RSGS and RMGS (Figs. 3c,g and 4g,i). Soil type (sand fraction) or drought sensitivity did not exert significant effects on drought recovery (Figs. 3g and 4i). To test whether these findings were spatially dependent, we additionally investigated the relationships between drought recovery and different bioclimatic factors across diverse combinations of vegetation types and Köppen–Geiger climate zones for both RSGS and RMGS. Consistent patterns in the relationships between drought recovery and different bioclimatic factors further confirmed that our findings were robust and not dependent on climate zones (Extended Data Figs. 5 and 6).

Discussion

The findings supported the hypothesis that vegetation phenology significantly impacted drought recovery. These impacts were strongly dependent on the interplay between vegetation phenology, drought timing and vegetation type. The effects of spring phenology on drought recovery included biophysical feedbacks, carbon and nutrient allocation and other seasonal biological processes11,12,15. We found that earlier spring phenology during the drought year postponed drought recovery when extreme drought events occurred in the mid-growing season (Extended Data Fig. 4). Earlier spring phenology probably led to a longer growing season and enhanced spring vegetation growth, which increased evapotranspiration and thus increased the drawdown of soil moisture, resulting in progressive water stress11,14,27. Consequently, this early-spring phenology-associated biophysical process probably caused more negative effects on the water supply for summer and autumn, a period when vegetation growth is most sensitive to water availability (Extended Data Fig. 7) and thus postponed vegetation recovery from mid- and late-growing season drought. This potential mechanism was supported by our findings that: (1) vegetation growth cannot fully recover from an extreme drought event within a single growing season for more than 60% of extreme drought events when they occur in the mid-growing season (Fig. 2); and (2) the spring phenology and drought recovery maintained a negative relationship when extreme droughts occurred in mid-growing season irrespective of vegetation types (Extended Data Fig. 4). Earlier spring phenology could also increase the risk of spring frost damage and further deplete the carbohydrate storage for woody plants, leading to a negative imbalance in carbon supply for sustaining drought recovery28. By contrast, we found that earlier spring phenology during the drought year shortened drought recovery, whereas delayed spring phenology lengthened drought recovery when extreme drought events occurred in the early-growing season (Extended Data Fig. 4). This was probably because the effects of spring phenology on vegetation physiology outweighed the biophysical effects on drought recovery. Earlier spring phenology during the drought year enhanced spring vegetation growth and led to stronger positive effects on subsequent vegetation growth15. The effects of phenology on soil moisture were probably weaker when extreme drought events occur in the early-growing season relative to the mid-growing season, due to weaker progressive water consumption from the start of spring to the beginning of the extreme drought event.

For drought recovery in the subsequent growing season (RMGS), delayed spring phenology in the subsequent year slowed down vegetation growth and led to longer drought recovery (Fig. 4a,i). This effect of spring phenology in the subsequent year appeared to be a major process governing drought recovery under RMGS for all vegetation types regardless of drought timing (Fig. 4 and Extended Data Figs. 4 and 8). We further found that climatic factors during the drought year only had a marginal effect on spring phenology in the subsequent year (Supplementary Text and Extended Data Fig. 8). The direct effects of climate anomalies in the drought year on RMGS were much less important than the effect of spring phenology in the subsequent year (Fig. 4i and Extended Data Figs. 4 and 8). Higher SWE and lower temperature in the non-growing season dormancy period further delayed drought recovery in the subsequent year (Fig. 4 and Extended Data Fig. 8). This finding was consistent with previous studies indicating that more snow cover and lower winter temperature were expected to delay snowmelt and spring phenology and thus vegetation growth, particularly at high latitudes20,29,30,31. Alternatively, less winter snow cover may decrease soil moisture availability and limit vegetation growth, particularly in the mid-latitudes32,33,34. In the mid-latitudes, climate warming had induced earlier snowmelt leading to drier soils in the growing season35, thereby decreasing vegetation growth and productivity. Additionally, low winter temperature would probably expose soil to freezing, which can prevent infiltration of snowmelt water and reduce spring and summer soil moisture supply. This reduction in spring and summer soil moisture in turn reduced vegetation growth and postponed subsequent growing season drought recovery. Low winter temperature could also increase frost injury thus reducing vegetation growth in the subsequent growing season36.

There was an important role of preceding growth condition on drought recovery (Figs. 3c and 4g). As expected, this effect was more prominent in woody plants (forests and shrubs), particularly when drought recovery extended into the subsequent growing season (Fig. 4g). On the one hand, more favourable preceding growth conditions were normally associated with higher drought resistance37,38. On the other hand, favourable preceding growth conditions could increase plant photosynthesis and produce more non-structural carbohydrates, which provided a vital carbon resource supporting more rapid vegetation postdrought recovery39,40,41. Better preceding growth condition could also stimulate root growth which can aid vegetation access to deeper groundwater resources and shorten drought recovery42. Divergent effects of preceding growth conditions on drought recovery among different vegetation types (Fig. 4g) were probably associated with the difference in anatomical structures of plant organs and the seasonal dynamics of non-structural carbohydrates. Nevertheless, the physiological basis of these divergent effects remained poorly understood.

The analyses further revealed that vegetation phenology effects were compounded by high temperatures. More specifically, increased mean annual temperature was significantly associated with a slowdown of drought recovery (Figs. 3d and 4h). This was probably due to the negative impacts of higher temperatures on plant water availability, which can lead to significant reductions in vegetation growth43. This finding implied that a projected future warmer climate could intensify the effects of spring phenology on vegetation recovery from extreme droughts, which were increasing in both frequency and severity21.

Although some studies indicated that the interannual climate variability may reshape the climate sensitivity of vegetation growth44,45, we did not observe significant effects of interannual variability of temperature and precipitation on drought recovery. We also did not find that drought sensitivity or drought response lags impacted drought recovery (Methods), although vegetation with higher drought sensitivity and shorter drought response lags was mainly distributed in warm dry regions (Extended Data Figs. 1 and 9) that were more subject to reduced levels of soil water availability and more vulnerable to drought6,46.

In summary, we reported two different trajectories in drought recovery, one occurring within a single growing season and another where recovery persists to the subsequent growing season. It showed that vegetation growth could not fully recover within a single growing season for more than 50% of extreme drought events. These findings highlighted crucial vegetation phenology effects on drought recovery, particularly emphasizing the important ramifications of shifts in spring phenology on drought recovery in the subsequent year. We highlighted that the importance of vegetation phenology effects on drought recovery was comparable to, or even larger than, other more well-known pre- and postdrought climatic factors. These findings added new insights toward better understanding terrestrial ecosystem resistance and resilience to extreme climate events. We further highlighted that the interactions between vegetation phenology and drought needed to be well represented in Earth system models to better predict the future trajectories of terrestrial ecosystems in a more extreme climate regime.

Methods

Vegetation growth proxies

Three complimentary vegetation growth proxies were used to quantify drought recovery, including the NDVI as a measure of vegetation greenness, CSIF as a metric of photosynthesis and microwave-based VOD as a measure of canopy biomass change. The third-generation NDVI, derived from the advanced very high-resolution radiometer (AVHRR), is provided by the Global Inventory Modelling and Mapping Studies group (GIMMS NDVI3g; https://climatedataguide.ucar.edu/climate-data/ndvi-normalized-difference-vegetation-index-3rd-generation-nasagfsc-gimms). The biweekly GIMMS NDVI3g has a spatial resolution of 0.083° and covers 1982–2015. For GIMMS NDVI3g, errors and noise effects have been corrected for inconsistencies resulting from atmospheric aerosols, sensor degradation, orbit drift, cloud contaminations, solar zenith angle, viewing angle effects due to satellite drift and volcanic aerosols22,47. GIMMS NDVI3g data have been widely used in monitoring changes in vegetation activity48,49, vegetation phenology50,51 and productivity52 both regionally and globally. Regions with multiyear mean annual NDVI values <0.2 during 1982–2015 were discarded from the final analyses. The global gridded CSIF dataset was obtained from https://osf.io/8xqy6/ (ref. 53). The CSIF dataset was trained by a machine learning approach based on surface reflectance from the moderate resolution imaging spectroradiometer (MODIS) aiming at reproducing SIF from the Orbiting Carbon Observatory-2. Two products are available, including clear-sky instantaneous SIF and the all-sky daily SIF. In this study, we used the all-sky daily SIF product, which covers 2000–2016 and has a 0.05° spatial resolution and 4 d temporal resolution. The microwave-based VOD data are an indicator of total vegetation water content in aboveground biomass and thus strongly connect to aboveground biomass, especially at seasonal timescales25. VOD is derived using microwave remote sensing and is thus mostly independent from the vegetation proxies based on visible and near-infrared remote sensing. The VOD is less affected by sun illumination and atmospheric effects than are optical vegetation indices, such as NDVI. The VOD is also sensitive to vegetation phenology54 and has been used to investigate soil moisture and drought-related impacts on vegetation growth25,55,56. In this study, the X-band (10.7 GHz) VOD was obtained from the advanced microwave scanning radiometer (AMSR) land parameter data record (LPDR, v.2, http://files.ntsg.umt.edu/data/LPDR_v2/)56. The LPDR was generated using calibrated, multifrequency daily brightness temperature retrievals from the AMSR for end of growing season (EOS) (AMSR-E, 2003–2010) and AMSR2 (2013–2015). The resulting LPDR dataset provides a long-term (June 2002–December 2015) gridded VOD record with a daily temporal resolution and a spatial resolution of 0.25°. In this study, we analysed VOD data during 2003–2015. The three vegetation growth proxies were resampled into a spatial resolution of 0.5° to match the other datasets.

Climate, drought and soil datasets

Gridded monthly climate data for the period 1982–2015 were obtained from the Climatic Research Unit (CRU TS4.02, http://www.cru.uea.ac.uk/cru/data/hrg/), including precipitation (Precip), temperature (Temp) and potential evapotranspiration (PET) at a spatial resolution of 0.5° (ref. 57). The vapour pressure deficit (VPD) was calculated based on the CRU climate datasets, following the method of Wu et al.7. The gridded monthly SWE during 1982–2015 with a spatial resolution of 25 km was obtained from the European Space Agency (ESA) GlobSnow dataset (https://www.globsnow.info/). GlobSnow combines SWE retrieved from multisatellite microwave observations with forward snow emission model simulations and ground-based weather station observations for non-mountainous regions of the Northern Hemisphere. Because of the improved accuracy achieved by assimilating independent sources of information, this is considered the best SWE product currently available for climate analysis58.

The SPEI59 was used as a drought metric in this study. SPEI is a simple and physiologically relevant drought index based on a water budget (difference between potential evapotranspiration and precipitation), which makes it more relevant to ecosystem water stress than meteorological drought indices based only on precipitation and temperature26. Besides, its multiscalar characteristics representing temporally lagged (for example, 1–24 month) integration times (for example, 3 month SPEI integrates water status over the previous 3 months) enable identification of different drought types and impacts48. Here, gridded monthly SPEI in the period 1982–2015 with a spatial resolution of 0.5° was extracted from the Global SPEI dataset (https://spei.csic.es). In our main analyses, we defined extreme droughts using SPEI at timescale of 3 months (SPEI3) because vegetation growth in the Northern Hemisphere responds predominantly to SPEI within timescales of 2–4 months7,48. The findings derived from different definitions of extreme droughts using SPEI at timescales of 1, 6, 9 and 12 months were further compared. Gridded soil sand content data were obtained from the Harmonized World Soil Database v.1.2 (https://www.fao.org/soils-portal/data-hub/soil-maps-and-databases/harmonized-world-soil-database-v12/en/)60. This database has 30 arcsec spatial resolution and >15,000 different soil mapping units that combine existing regional and national updates of soil information worldwide.

Climate zones and vegetation types

Six major climate zones were defined from the Köppen–Geiger climate classification61, including arid (AR), warm dry (WD), warm humid (WH), cold dry (CD), cold humid (CH) and polar regions (Supplementary Table 1 and Supplementary Fig. 4). The classifications of four vegetation types, including evergreen forests, deciduous forests, shrubs and grasses, are based on the GLC 2000 land use/land cover dataset from the European Commission Joint Research Center (https://forobs.jrc.ec.europa.eu/products/glc2000/glc2000.php). The land cover map discriminates 17 land cover types on the basis of the IGBP (International Geosphere–Biosphere Programme) classification scheme. In this study, the evergreen broadleaf forests and evergreen needleleaf forests were grouped into a single ‘evergreen forests’ class; and the deciduous needleleaf forests and deciduous broadleaf forests were grouped into a ‘deciduous forests’ class. The categories of closed shrublands, open shrublands, woody savannas and savannas were grouped into a single ‘shrubs’ class. Grasslands are denoted as ‘grasses’. The regrouped vegetation map was further resampled to a 0.5° spatial resolution to match the other datasets (Supplementary Fig. 4).

Vegetation phenology

The HANTS-Maximum method was applied to filter and interpolate the biweekly NDVI into daily time series62. Then, a threshold-based method was used to retrieve the start of growing season (SOS) and the end of growing season (EOS). This method has been widely used and proven to be an effective tool to capture the vegetation phenology at large scale63,64. In this study, the analyses were concentrated on regions with a single growing season occurring in each calendar year, featured by a nearly normal NDVI distribution with a summer peak. Therefore, we first discarded pixels having multiple growing seasons occurring in a single year or without clear phenological cycles. Then, we extracted the SOS and EOS dates for each pixel using a pixel-specific threshold of minimum value plus 30% of the seasonal amplitude fitting for multiyear averaged smoothed NDVI24 (Supplementary Fig. 5). The growing season was then defined as the period from SOS to EOS for each pixel. The mid-growing season (mid-GS) was defined as the two consecutive months with the largest NDVI value but occurring no earlier than April or later than October15. The early-growing season (early-GS) was defined as the period from the month when SOS occurred to the beginning of the mid-GS and the late-growing season (late-GS) is defined as the period from the end of the mid-GS to the EOS. The dormant season was defined as the period from the EOS of the current year to the SOS of the next year. Overall, about 81% (NDVI-based) of vegetated land was analysed following this screening process.

Identification of drought events

An extreme drought event is defined as a period with SPEI3 < −2 and reduced vegetation greenness and productivity (lower average NDVI during drought period than the long-term mean value in the period 1982–2015) in early-GS, mid-GS and late-GS, respectively. An extreme drought event ends when SPEI3 > −2 (ref. 4). Only single extreme drought events within one growing season are considered in this study to avoid longer legacy effects on recovery from repeated extreme climate events65; another extreme climate event (either extreme drought or extreme wetness) occurring within 4 yr before and after a given drought event was not included in this analysis. Extreme drought events were alternatively defined by the SPEI at timescales of 1, 6, 9 and 12 months. The drought recovery from different extreme drought definitions using SPEI at a variety of timescales was then quantified and compared.

Drought response and drought recovery

The drought response was quantified by two parameters: drought response lag and drought sensitivity. Drought response lag was measured by the duration (months) between the minimal SPEI3 anomalies and the deepest suppression of NDVI induced from the drought events. Drought sensitivity was determined by the Pearson’s correlation coefficient between mean growing season NDVI and SPEI3 in the period 1982–2015. Both drought response lag and drought sensitivity were included into the RF model to quantify their effects on drought recovery. Pearson’s correlation analysis was additionally performed to evaluate the responses of seasonal mean NDVI to SPEI3 in early-, mid- and late-GS. The seasonal mean NDVI and SPEI3 were linearly detrended before the Pearson’s correlation analysis.

Drought recovery is defined as the duration (months) starting from the month with the deepest suppression of NDVI to the month when NDVI returns to within 95% of the long-term average baseline. The monthly SPEI3 and NDVI time series were first smoothed by a 3 month forward moving window. They were then sequentially deseasonalized and linearly detrended. The multiyear mean value of the detrended NDVI in the period 1982–2015 was calculated as the long-term average baseline for each pixel. Thus, at each pixel, there is a location-specific, stationary baseline for the whole period 1982–2015. In this study, two different drought recovery trajectories were identified, including: (1) vegetation recovering in the same year of a drought event and before the dormant season (RSGS); and (2) vegetation recovery extending through the dormant season and into subsequent year (RMGS; Fig. 1a and Supplementary Fig. 1). To avoid lengthening the drought recovery duration due to algorithm design, in the RMGS, the drought recovery was calculated as the total length of the recovery period minus the length of the dormant season. Severe drought events and subsequent drought recovery were pixel specific and calculated for every combination of NDVI and SPEI3. Only drought and recovery events that were fully contained in the period 1982–2015 were analysed. The drought recovery derived from CSIF and VOD using the same protocol was additionally compared (Extended Data Fig. 3). The comparison of three different vegetation metrics allows for validating the NDVI-based findings and confirming whether the drought recovery patterns are consistent for different vegetation proxies.

Factors affecting drought recovery

Separate RF models were built for attributing drought recovery to vegetation phenology, different bioclimatic and soil factors under either RSGS (with sample size of 29,614 events) or RMGS (with sample size of 42,733 events). For drought recovery under RSGS, 15 variables were included as predictive factors. These factors can be stratified into four groups representing unique processes and properties. First, the start of the growing season in the drought year (SOSdrought_year) was introduced into the model to capture the spring phenology effects on drought recovery. Second, multiple climatic variables during the predrought period, drought response lag period and postdrought recovery period were considered, including the precipitation (Precippre, Preciplag and Precippost, respectively) and mean VPD (VPDpre, VPDlag and VPDpost, respectively). The mean NDVI anomaly in the 6 months within growing season before an extreme drought event (NDVIpre) was also included in the model. Moreover, the effects of long-term climate conditions were considered, including mean annual precipitation (MAP), mean annual temperature (MAT) and the multiyear mean water balance (calculated as mean annual precipitation minus mean annual potential evapotranspiration) and the variance of hydroclimatic conditions, including the coefficient of variation (CV) of annual precipitation (PrecipCV) and annual temperature (TempCV) in the period 1982–2015. Third, drought sensitivity, which capture the response of vegetation growth to drought, were also included into the final model. Fourth, the soil sand content was incorporated. For the drought recovery under RMGS, we additionally considered the effects of SOS in the subsequent year (SOSsubsequent_year), the length of dormant period (dormant length), the SWE and the mean temperature (Tempdormant) in the dormant period and drought response lag. Anomalies in these bioclimatic variables in different periods corresponding to individual extreme drought event were calculated on the basis of the deseasonalized and linearly detrended series.

Three complementary approaches were used to determine which factors were most important for drought recovery. First, the Pearson’s correlation between explanatory variables and the drought recovery was calculated. Second, to avoid multicollinearity among the variables, a matrix of pairwise correlation was calculated and any variable with high correlation (R > 0.5) with other predictor variables was removed. Each pairwise correlation was performed and the variable with the lower correlation with the dependent variable was removed. Finally, the RF model was applied to estimate the importance of each explanatory variable.

The RF extends the standard classification and regression tree (CART) method by creating a collection of classification trees through binary divisions66,67. Here, we defined that each RF had 500 binary trees with one covariate chosen at random from the full set to determine the splitting rule, with a minimal terminal node size of 5. With the RF mode, variable importance ranking was also calculated by the permutation importance method. Permutation feature importance is a model inspection technique used for any fitted models, which is especially useful for nonlinear models. The importance of permutation features is defined as the decrease in the model accuracy when a single feature value is randomly scrambled. The larger the decrease after the permutation, the greater the importance of this variable in the forest and in the prediction accuracy68. The response function of drought recovery to each individual factor independent of other covariates is shown as a partial dependent plot. Partial dependence gives the marginal effect of a covariate on the response variable. The bootstrap method was used to compute the average over all the predictions of drought recovery and the 95% confidence interval (CI).