Widespread spring phenology effects on drought recovery of Northern Hemisphere ecosystems

The time required for an ecosystem to recover from severe drought is a key component of ecological resilience. The phenology effects on drought recovery are, however, poorly understood. These effects centre on how phenology variations impact biophysical feedbacks, vegetation growth and, ultimately, recovery itself. Using multiple remotely sensed datasets, we found that more than half of ecosystems in mid- and high-latitudinal Northern Hemisphere failed to recover from extreme droughts within a single growing season. Earlier spring phenology in the drought year slowed drought recovery when extreme droughts occurred in mid-growing season. Delayed spring phenology in the subsequent year slowed drought recovery for all vegetation types (with importance of spring phenology ranging from 46% to 58%). The phenology effects on drought recovery were comparable to or larger than other well-known postdrought climatic factors. These results strongly suggest that the interactions between vegetation phenology and drought must be incorporated into Earth system models to accurately quantify ecosystem resilience. The authors reveal complex drought recovery responses to phenology shifts, in that early spring can shorten or lengthen recovery, while delayed spring following drought events delays it. These effects suggest a need to incorporate phenology aspects into resilience models.

Article https://doi.org/10.1038/s41558-022-01584-2 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 drought 21 . Changes in vegetation greenness caused by droughts are identified by satellite observations of normalized difference vegetation index (NDVI;1982 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 25 (Methods). Extreme drought events are defined using the multiscalar standardized precipitation evapotranspiration index (SPEI 26 ; 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: R SGS and R MGS ( 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 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 seasons 14 . 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 drought 15 .
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 cycle 12,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 (R SGS ) or persist to subsequent growing seasons (within multiple growing seasons) (R MGS ) (Methods). Vegetation phenology interacts with drought timing and thus determines subsequent vegetation growth 11,14,15 (Fig. 1). On the one hand, spring phenology strongly governs vegetation growth and soil water consumption in the growing season [11][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 processes 19,20 , thus influencing drought recovery (Fig. 1). Vegetation phenology is already being rapidly impacted by climate warming 12,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). 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 R SGS and R MGS 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). Overall, drought recovery following early-growing season droughts was on average longer than that following mid-and late-growing season droughts. Under R SGS , 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 R MGS , 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 R SGS and R MGS (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 R SGS and R MGS . 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 (R SGS and R MGS ) and four vegetation types (evergreen forests, deciduous forests, shrubs and grasses). The RF models performed well in capturing drought recovery, with better performance for R MGS (R 2 ranged from 0.86 to 0.90 across vegetation types) than for R SGS (R 2 ranged from 0.71 to 0.75).
The analyses identified significant but differentiated phenological effects on drought recovery among different vegetation types and between R SGS and R MGS . For R SGS , 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 R SGS , 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).  We next tested the second hypothesis that a delayed spring phenology in the subsequent year will postpone drought recovery. For R MGS , 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 R MGS 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 R SGS and R MGS (Figs. 3a,b and 4b,c). Higher mean annual temperature led to longer drought recovery for both R SGS and R MGS , 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 R SGS or R MGS . The interannual variability of mean annual temperature and mean annual precipitation also showed no significant effects on drought recovery for either R SGS or R MGS (Figs. 3g and 4i).
Finally, we found that the NDVI pre (mean NDVI anomaly in the 6 months within growing season preceding extreme drought) exerted important effects on drought recovery, with more prominent effects in R SGS (with importance of 25-31% across vegetation types). For instance, negative anomalies in NDVI pre tended to delay drought recovery under both R SGS and R MGS (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 R SGS and R MGS . 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 processes 11,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 stress 11,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 recovery 28 . 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 year , d) (e) and mean VPD during the drought response lag period (VPD lag ) (f). The y axis in a-f shows the drought recovery (months). Shaded areas of different colours in a-f 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). Precip pre , total precipitation in the 6 months preceding the extreme drought events; VPD pre , mean VPD in the 6 months preceding the extreme drought events; Precip lag , total precipitation during the drought response lag period; MAP, mean annual precipitation; Precip CV , coefficient of variation of mean annual precipitation; and Temp CV , coefficient of variation of mean annual temperature. Numbers in the legend indicate the sampling sizes. The z-scores were calculated for Precip post , VPD post , NDVI pre , VPD lag , Precip pre , VPD pre and Precip lag before the RF analysis.
Article https://doi.org/10.1038/s41558-022-01584-2 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 growth 15 . 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 (R MGS ), 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 R MGS 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 R MGS 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 latitudes 20,29-31 . Alternatively, less winter snow cover may decrease soil moisture availability and limit vegetation growth, particularly in the mid-latitudes 32-34 . In the mid-latitudes, climate warming had induced earlier snowmelt leading to drier soils in the growing season 35 , thereby decreasing vegetation growth and productivity. Additionally, low winter temperature would probably expose soil to freezing, which can Shaded areas of different colours in a-h 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 a-h 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). Precip pre , total precipitation in 6 months preceding the extreme drought events; SOS drought_year , the start of growing season of drought year; MAP, mean annual precipitation; Precip CV , coefficient of variation of mean annual precipitation; Temp CV , coefficient of variation of mean annual temperature; VPD pre , mean VPD in 6 months preceding the extreme drought events; Precip lag , total precipitation during the drought response lag period; and VPD lag , mean VPD during the drought response lag period. Numbers in the legend indicate the sampling size. z-scores were calculated for Precip post , VPD post , SWE dorm , Temp dorm , NDVI pre , Precip pre , VPD pre , Precip lag and VPD lag before the RF analysis.
Article https://doi.org/10.1038/s41558-022-01584-2 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 season 36 . 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 resistance 37,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 recovery 39-41 . Better preceding growth condition could also stimulate root growth which can aid vegetation access to deeper groundwater resources and shorten drought recovery 42 . 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 growth 43 . 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 severity 21 .
Although some studies indicated that the interannual climate variability may reshape the climate sensitivity of vegetation growth 44,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 drought 6,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.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41558-022-01584-2. Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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 aerosols 22,47 . GIMMS NDVI3g data have been widely used in monitoring changes in vegetation activity 48 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 timescales 25 . 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 phenology 54 and has been used to investigate soil moisture and drought-related impacts on vegetation growth 25, 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  Because of the improved accuracy achieved by assimilating independent sources of information, this is considered the best SWE product currently available for climate analysis 58 .
The SPEI 59 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 temperature 26 . 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 impacts 48 . 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 months 7,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 classification 61 , 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 series 62 . 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 scale 63,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 NDVI 24 (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 October 15 . 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 Nature Climate Change Article https://doi.org/10.1038/s41558-022-01584-2 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 events 65 ; 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 (R SGS ); and (2) vegetation recovery extending through the dormant season and into subsequent year (R MGS ; Fig. 1a and Supplementary Fig. 1). To avoid lengthening the drought recovery duration due to algorithm design, in the R MGS , 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 R SGS (with sample size of 29,614 events) or R MGS (with sample size of 42,733 events). For drought recovery under R SGS , 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 (SOS drought_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 (Precip pre , Precip lag and Precip post , respectively) and mean VPD (VPD pre , VPD lag and VPD post , respectively). The mean NDVI anomaly in the 6 months within growing season before an extreme drought event (NDVI pre ) 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 (Precip CV ) and annual temperature (Temp CV ) 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 R MGS , we additionally considered the effects of SOS in the subsequent year (SOS subsequent_year ), the length of dormant period (dormant length), the SWE and the mean temperature (Temp dormant ) 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 divisions 66,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 accuracy 68 . 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).

Code availability
All codes 69 used to analyse data and create figures are available through the GitHub at https://github.com/leeyang1991/ phenology-effects-on-drought-recovery.  Fig. 2 | Drought recovery as a function of aridity index over the mid-and high-latitudinal Northern Hemisphere. Aridity index (AI) is calculated as the ratio of annual total precipitation to annual total potential evapotranspiration using the gridded climate data of CRU TS4.02 (http://www. cru.uea.ac.uk/cru/data/hrg/). The spatial pattern in the multiyear mean aridity index in the period of 1982-2015 was shown in a. The relationship between drought recovery (averaging of both of that within a single growing season and that within the subsequent growing season) and multiyear mean aridity index in the mid-and high-latitudinal Northern Hemisphere was shown in b. Blue bars in b denote the percentage of area with corresponding AI interval to the total study region. Shaded area in b indicates the 95% confidence interval for the average drought recovery in different AI bins. According to the classification of UNEP (1992), areas with AI > 0.6, 0.65-0.5, 0.5-0.2, 0.2-0.05, and < 0.05 are defined as humid region, sub-humid region, semi-arid region, arid region, and hyper-arid region, respectively. Contribution of different factors was estimated by random forest models dependent on drought timing and vegetation types for drought recovery occurring within a single growing season (a, R SGS ) or within the subsequent growing season (b, R MGS ), respectively. Red and blue colours in a and b indicate positive and negative correlation coefficients between predicting variables and drought recovery, and darker colours indicate stronger correlation. Only factors with contribution >2% (averaged contribution of all factors) are displayed. The predrought period represents the six months preceding the extreme drought events. We did not consider the vegetation recovery under late-growing season drought for the R SGS scenario, to prevent the problem of overfitting caused by insufficient training dataset samples. We regrouped the extreme drought events into early-growing season (Early-GS) droughts, mid-growing season (Mid-GS) droughts, and late-growing season (Late-GS) droughts. Four kinds of vegetation types were considered, including evergreen forests (Evergreen), deciduous forests (Deciduous), shrubs, and grasses. A variety of bioclimatic and edaphic factors were introduced into this analysis, including start of growing season of drought year (SOS drought_year ), start of growing season of subsequent year (SOS subsequent_year ), total precipitation in preceding drought period (Precip pre ), postdrought period (Precip post ) and drought response lag period (Precip lag ), vapour pressure deficit in preceding drought period (VPD pre ), postdrought period (VPD post ), and drought response lag period (VPD lag ), NDVI during preceding drought period (NDVI pre ), dormant length, mean temperature and snow water equivalent during dormancy period (Temp dorm and SWE dorm ), both mean values and variability of mean annual temperature (MAT and Temp CV ) and mean annual precipitation (MAP and Precip CV ), multiyear mean water balance, drought sensitivity, drought response lag, as well as sand fraction. Four different kinds of vegetation types were considered, including evergreen forests (Evergreen), deciduous forests (Deciduous), shrubs and grasses. In this analysis, we divided our study region into six major Köppen-Geiger climate zones, including arid (AR), warm humid (WH), warm dry (WD), cold humid (CH), cold dry (CD) and polar regions (Table S1). Lines in a-f are the linear fits. Error bars for each data point in a-f represent the standard deviations of the data in different combinations of climate zones and vegetation types. Grey shaded areas in a-f represent the 95% confidence intervals for the linear fits. The sample number is 24 for each linear fit. Six important factors were considered in this analysis, including total precipitation in the postdrought recovery period (Precip post ), the mean vapour pressure deficit in both postdrought recovery period (VPD post ) and drought response lag period (VPD lag ), the mean NDVI in six months preceding the extreme drought (NDVI pre ), mean annual temperature (MAT, °C) and the start of growing season in drought year (SOS drought_year , days). The direct effects of climate anomalies during drought year (including mean temperature, Temp drought_year , total precipitation, Precip drought_year , and mean vapour pressure deficit, VPD drought_year ), during the dormant period (including the mean temperature, Temp dormant , and snow water equivalent, SWE dormant ), and during the spring of the subsequent year (including mean temperature, Temp subsequent_year_spring , and total precipitation, Precip subsequent_year_spring ), on drought recovery within the subsequent growing season were considered. The indirect effects of Temp drought_year , Precip drought_year , VPD drought_year , Temp dormant , and SWE dormant on the spring phenology of subsequent year and then the drought recovery was considered. The arrows indicate the hypothesized direction of causation, with positive and negative relationships in blue and red, respectively. Coloured lines represent significant (p < 0.05) relationships by the two-sided t-test. Line thickness is proportional to the strength of the relationship and to the standardized path coefficients adjacent to each arrow. We used the χ 2 test (and associated p values), goodness-of-fit index (GFI), an adjusted goodness-of-fit index (AGFI), and the root mean square error of approximation (RMSEA), to evaluate the fit of the established structural equation model. The SEM model is with a χ 2 = 0.002, RMSEA = 0.000, GFI=1.000, and AGFI=1.000.