Global Assessment of Cumulative and Time-Lag Effects of Drought on Land Surface Phenology

ABSTRACT Increased frequency and intensity of droughts under climate change will have a significant impact on land surface phenology, however, the drought-phenology interactions that are associated with complex temporal effects are not well understood. This study examined the response of land surface phenology to drought cumulative and time-lag effects by using the standardized precipitation and evapotranspiration index (SPEI) and explored the influence of hydrothermal and plant physiological factors on the phenology-drought relationship. We used the maximum correlation (rmax-cml ) between phenology and cumulative SPEI (1- to 12-month) to determine the cumulative effect. The maximum correlation (rmax-lag ) between phenology and lagged SPEI (1-month) was utilized to analyze the time-lag effect. Overall, the cumulative effect affected 25.36% of the vegetated area at the start of the growing season (SOS), 26.43% of the end of the growing season (EOS), and 26.57% of the length of the growing season (LOS). SOS was negatively affected by long-term SPEI, whereas EOS and LOS had positive correlations with short-term SPEI. The rmax-cml for shrubland was the largest, and the SOS and LOS response time scales of the forest were the shortest. The rmax-cml was larger in arid and semi-arid (AR and SAR) than in humid and semi-humid (HU and SHU). Meanwhile, the response time scales were longer in HU and SHU than in AR and SAR. The time-lag effect had a larger area of impact on land surface phenology than the cumulative effect, the areas were 46.12%, 47.93%, and 50.45% for SOS, EOS, and LOS, respectively, and the lagged time scales were longer. The phenology-SPEI correlation was dominantly driven by hydrological conditions, and the time scales were mainly affected by thermal factors. Moreover, the onset of phenology and the growth/senescence rate of plants influenced the relationship, suggesting that hydrothermal conditions may shift the time of phenology by regulating the growth/senescence rate and may thus modulate the phenology–drought interaction.


Introduction
Global climate change has disrupted many ecosystem processes and functions in recent decades, including carbon and water cycles and energy flow (Zheng et al. 2018), and dramatically affected vegetation productivity, biodiversity, and land surface phenology (Chen et al. 2014;Zhang et al. 2018b;An et al. 2020). The timing of land surface phenology events is an important signal for the inception and cessation of carbon uptake during the ecosystem carbon cycle (Barichivich et al. 2012). Given the intense and observable connection between plant life cycle events and climate variables in terrestrial ecosystems, land surface phenology has elicited much attention. Previous research, however, has primarily concentrated on the start of the growing season (SOS), with a few studies focusing on the end of the growing season (EOS) and the length of the growing season (LOS) (Uttaruk and Laosuwan 2017;Dannenberg et al. 2015). Growing evidence suggests that vegetation phenology has an important contribution to the control of ecosystem productivity, and water and carbon cycling (Zarei, Shabani, and Mahmoudi 2020;Penuelas, Rutishauser, and Filella 2009), as a result, understanding the feedback relationship between climate change and land surface phenology is important for gaining a better understanding of the critical processes that manage ecosystems.
The dynamics of land surface phenology and associated environmental variables could detail the spatiotemporal variability of vegetation during the terrestrial biogeographic cycle, and mitigation strategies for climate change (Cui and Shi 2021;Li et al. 2018;Piao et al. 2015). It has been demonstrated in numerous studies that changes in the phenology of surface plants are influenced by hydrothermal conditions such as temperature, precipitation, and radiation Zu et al. 2018). In temperate grasslands, the growing season was dominated by temperature and precipitation, and early precipitation has a remarkable effect on autumn phenology Liu et al. 2016). In regions with mid-tohigh latitudes, temperature as the main factor of land surface phenology has long been recognized, and growing season indicators based on temperature (or heat) have been introduced to assess the indirect effects of climatic warming on land biomes (Cui et al. 2018;Barichivich et al. 2013). In addition, a growing body of research suggests that global warming has induced early-start, late-end, and prolonged growing seasons in recent decades (Ren et al. 2019). The complex interaction between plants and the environment is apparent in these studies concerning how land surface phenology responds to environmental changes. Consequently, further insightful and systematic analyses are needed to enhance the understanding of these interactions.
Drought is a more complex abiotic stress affecting terrestrial ecosystems (Vicente-Serrano et al. 2012;Laosuwan et al. 2016). Drought stresses accompanied by higher frequencies and intensities are expected to pose significant risks to the biodiversity of terrestrial ecosystems as well as to the stability of producing agriculture Zhou et al. 2019;Zarei and Mahmoudi 2020). Therefore, the rapid and precise estimation of drought stressinduced declines in vegetation activity is a great challenge in terrestrial ecosystem research for the development of mitigation strategies for zones vulnerable to dryland disturbance. However, the mechanism of phenology-drought feedback is still controversial as the response of surface phenology to drought stress is associated with intricate temporal effects (Yun et al. 2018;He et al. 2018). Bernal, Estiarte, and Penuelas (2011) suggested that SOS exhibits an advanced trend in drought conditions, whereas Kang, Wang, and Liu (2018) found that drought may delay SOS. Moreover, drought advances EOS in the grasslands of China and Canada in summer and autumn (Cui, Martz, and Guo 2017), but it results in late EOS in the northeastern U.S for most deciduous forests (Xie, Wang, and Silander 2015;Xie et al. 2018). Furthermore, vegetation growth may be influenced by previous drought events rather than current moisture conditions. Plants often require time after a severe drought to rehabilitate their impaired root systems and recover their pre-drought development capabilities (Rammig et al. 2015). For effective management and conservation of terrestrial biomes, it is essential to investigate the response of land surface phenology variations to drought stress and mitigation strategies, as well as an important process to understand the intricate mechanisms of climatevegetation interactions in ecosystems.
Numerous studies have established that changing climate could affect the growth and dynamics of terrestrial plant communities (Nemani et al. 2003) through a complicated temporal effect (Anderegg et al. 2015). Specifically, there is an asymmetric way in which climate dynamically affects vegetation, i.e. with or without time-lag (Wu et al. 2015;Wen et al. 2018) and cumulative effects (Ivits et al. 2016;Zhang et al. 2018a). Therefore, exploring the cumulative and time-lag effects induced by water accumulation values at different periods is of importance for the comprehension of the complex drought-plant relationships and for implementing effective vegetation management strategies. The time-lag effect indicates the impact of previous drought stress occurring at certain time points on the present plant growth (Wu et al. 2015), and the cumulative effects of meteorological conditions suggest vegetation growth is heavily dynamically influenced by them (e.g. precipitation and evapotranspiration) of the previous and current months (Vicente-Serrano et al. 2013;Zhao et al. 2020). Given this relation, we can reasonably speculate that drought cumulative and time-lag effects have been essential in regulating land surface phenology. In addition, the temporal effect on vegetation from climate varies with vegetation type and climatic conditions, i.e. the same climatic conditions have distinct temporal effects on different biomes, despite having different temporal effects on the same biomes with different climatic conditions (Wu et al. 2015). It is difficult to predict how vegetation growth will respond to climate change using dynamic vegetation models given the limitations of understanding these temporal effects, thereby increasing the challenge of examining temporal variability in drought response.
Further effort should be devoted to examining climate-driven effects on vegetation's response to drought to effectively assess and predict vegetation dynamics.
Changes in environmental conditions and land biomes play a vital role in determining the difference in the responses of vegetation phenology to drought (Xu, Wang, and Zhao 2021;Kang, Wang, and Liu 2018). However, the impact of predominant hydrothermal conditions on the temporal effects of different plant life cycles at different growth stages under different land biomes and water gradients has not been scientifically investigated. This investigation is particularly valuable because of climatic and physiological factors during the plant growing season in terrestrial ecosystems. Our main objectives are to (a) quantify the cumulative and time-lag effects of recent droughts  on global land surface phenology, (b) investigate how these effects vary with land biomes and water gradients, and (c) analyze the effect of hydrothermal conditions and plant physiological factors on these effects.

Datasets
The land surface phenology products, including SOS, the rate of the greening season (ROG), EOS, the rate of the senescence season (ROS), and LOS, were obtained from the Vegetation Index and Phenology (VIP) Laboratory of the University of Arizona (https://vip. arizona.edu/viplab_data_explorer.php). The SPEI index was from the SPEIbase v.2.5 at the Consejo Superior de Investigaciones Científicas (CSIC, http:// digital.csic.es/handle/10261/153475), was used to determine the cumulative and time-lag effects on land surface phenology. The IGBP land cover classification dataset from the Land Processes Distributed Active Archive Center (LPDAAC, http://modis-land. gsfc.nasa.gov/landcover.html/) for classifying different vegetation types, in this study, was reclassified into 4 land biomes: Forest, Scrubland, Grassland and Cropland (Table S1 and Fig. S1). The Aridity Index (AI), calculated according to the United Nations Environment Programme (UNEP) criteria, was used to classify the different moisture gradients, including arid (AI < 0.20, AR), semi-arid (0.20 < AI < 0.50, SAR), sub-humid (0.50 < AI < 0.65, SHU) and humid (AI > 0.65, HU) zones (Fig. S2). Meteorological data were derived from the ERA5-Land dataset product (https:// cds.climate.copernicus.eu/cdsapp#!/search?type= dataset) of the European Center for Medium-range Weather Forecasts (ECMWF). The specific calculations and more details for the different datasets can be viewed in the Supplementary Material.

Data Analysis
The maximum correlation (r max-cml ) between different phenological parameters and the cumulative SPEI (1to 12-month) and the time scales corresponding to its occurrence was used to calculate the impact of the drought cumulative effect on land surface phenology. And the time-lag effect was determined by the maximum correlation (r max-lag ) between the land surface phenology and the lagged SPEI (1-month) and the corresponding time scales. Sen's analysis and Mann-Kendall (M-K) test were used to identify the spatiotemporal dynamics of surface phenology. To investigate the influencing factors of phenology-drought relationships, correlation analysis, and redundancy analysis were used. The data analysis was done in Matlab 2016, ArcGIS10.7, and Canoco 5.0 software, and the figures were plotted using Origin2021b. The specific analysis method formulas and details can be viewed in the Supplementary Materials, and the details of our framework are shown in Figure 1.

Responses of Land Surface Phenology to Drought at Cumulative Time Scales
The water balance in different periods simultaneously promotes or inhibits plant growth and development, which in turn affects the prevailing phases of plant growth and development (Zu et al. 2018;Liu et al. 2016). Therefore, the cumulative value of available water in each stage of the plant life cycle may have a cumulative effect on the next stage, thereby inducing a cumulative effect on the land surface phenology. The r max-cml between different phenological parameters and cumulative SPEI (1-to 12-month) and the percentages of positive and negative correlations on different cumulative time scales identified the cumulative effects on land surface phenology ( Figure 2). The cumulative months corresponding to the emergence of r max-cml between SOS and cumulative SPEI were mainly distributed in 12, 1, and 2 months, and the peak value was observed at the 12month (positive in 4.72% and negative in 8.35%) ( Figure 2a). The largest positive r max-cml was found at the 10-month (r = 0.428), while the absolute largest negative r max-cml was at the 4-month (r = −0.434) ( Figure 2b). Generally, the proportion of SOS showing a positive correlation with cumulative SPEI was smaller than the proportion of negative correlation at different cumulative months. The largest proportions of the cumulative months corresponding to the r max-cml occurrences of EOS and cumulative SPEI happened mainly in 1-4 months, and the maximum value was detected in the 1-month with a positive of 9.29% and a negative of 8.23% (Figure 2c). The largest positive r max-cml was found at the 1-month (r = 0.437), while the absolute largest negative r max-cml was at the 5-month (r = −0.434). Overall, the proportion of EOS showing a significant positive correlation (58.16%) with cumulative SPEI at different cumulative time scales was greater than the proportion of significant negative correlation (41.84%). The proportion of the cumulative months corresponding to the r max-cml between LOS and cumulative SPEI was mainly distributed in 11-, 12-and 1-month scales, accompanied by a maximum observed in 1-month (positive in 11.47% and negative in 4.32%) (Figure 2e). The largest positive r max-cml was observed at the 8-month (r = 0.445), while the absolute largest negative r max-cml was found at the 7-month (r = −0.434) ( Figure 2f). Generally, the percentage of positive correlation (60.30%) between LOS and SPEI was larger than that of the negative correlation (39.70%). The dominant time scales of phenology responses to drought differed due to the attributes of phenology. SOS was mainly affected by long-term drought, in contrast to EOS and LOS, which were primarily affected by short-term drought. Peng et al. (2019) examined the relationship between EOSdrought in the northern hemisphere from 1982 to 2015 and obtained a similar conclusion to ours, that is, EOS is mainly controlled by short-term drought.

Spatial Patterns of Drought Cumulative Effects on Land Surface Phenology
The spatial distribution of r max-cml between SOS and cumulative SPEI ( Figure 3a) and the cumulative months ( Figure 3b) corresponding to their occurrence revealed the influence of cumulative effects on SOS. The positive correlations (9.83%) were mostly distributed in southwestern Eurasia, Mongolia, northern Australia, central Africa, most of North America, and western South America. Negative correlations (15.53%) were primarily observed in most of northern Eurasia, central North and South America, and northern and southern Africa. The time scales where r max-cml occurred were mainly 12-, 1-, and 2-month, among which the 12-month had the maximum area (13.07%), followed by the 2-month (11.26%) and 1-month (11.00%) scales. Contrary to the SOS result, the proportion of vegetation cover areas showing a positive correlation between EOS and cumulative SPEI was higher than that of negatively correlated areas ( Figure 3c). The significant positive correlation areas accounted for 15.37% of the vegetated area and were mainly distributed in northern Russia, southwestern Europe, southern China, eastern Australia, northern and eastern Africa, most of southwestern North America, and eastern South America. Moreover, 11.06% of the vegetated area was negatively correlated and mostly distributed in central Asia and Europe, central Africa, central South America, and northwestern North America. The time scales corresponding to r max-cml were mainly concentrated in 1-4 months and accounted for 52.52% of the significant correlation area (Figure 3d). The 1-month scale occupied the largest area (17.51%). Meanwhile, approximately 12.77% and 12.30% of the maximum correlations were occupied by the 2-and 3-month SPEI, respectively. Overall, 16.02% of the vegetation cover areas showed a positive correlation between LOS and cumulative SPEI, which were mainly concentrated in northern Russia, southwest China, eastern Australia, most regions of Africa and South America, and most regions of western and southeastern North America. Significantly negatively correlated regions accounted for 10.55% and were located in central Eurasia, Western Australia, and northwestern North America ( Figure 3e). The time scales corresponding to the occurrence of r max-cml had large variability in space. Specifically, 15.78% of the vegetation area was affected by the 1-month SPEI, followed by 12-month (12.76%) and 11-month (10.50%) SPEI (Figure 3f). The interannual and spatial dynamics of land surface phenology from 1982 to 2015 indicated that SOS was trending advanced, EOS was delayed and LOS was prolonged (Figs. S3 and S4), and previous studies also demonstrated this trend change in land surface phenology (Bernal, Estiarte, and Penuelas 2011;Xie et al. 2018;Chen et al. 2020). There were diverse correlations between different phenology indicators and cumulative SPEI in this study, suggesting that cumulative water availability plays an important role in vegetation growth and development in terrestrial ecosystems.

Drought Cumulative Effect on Land Surface Phenology across Different Land Biomes and Water Gradients
The r max-cml and cumulative months showed different distribution patterns between land surface phenology and 1-to 12-month SPEI across land biomes ( Figure 4).
The r max-cml between phenology and SPEI indicates the degree to which phenology is affected by drought, and the corresponding cumulative months indicate the sensitivity of phenology to drought (i.e. the length of time for phenology to respond to drought) (Xu et al. 2018;Liu et al. 2021). The proportion of significant correlations among land biomes showed that forest had the smallest proportion, where all phenology types responded to the cumulative effect, and the proportion of shrubland and grassland land surface phenology significantly affected by the cumulative effect was higher than that of forest and cropland (Figure 4), which illustrates the important regulation of accumulated available water in mitigating the response of terrestrial ecosystems to drought stress. Moreover, the trend of r max-cml between all phenology and cumulative SPEI was consistent across land biomes, with the highest correlation for shrubland and the lowest for forest and cropland. It showed that the cumulative effect on the land surface phenology of shrubland and grassland was stronger than that of forest and cropland. The corresponding cumulative months varied considerably among the different land biomes. The response time scales of shrubland in SOS was the longest (6.63 months) and that of the forest was the shortest (5.93 months) (Figure 4b). The cumulative months of the forest in EOS was the longest (5.39 months) and cropland was the shortest (4.99 months) (Figure 4d). The grassland had the longest response cumulative months (6.71 months), and forest had the shortest (5.98 months) for LOS (Figure 4f).
The fundamental mechanism is probably connected to the timing of land surface phenology and the biochemical characteristics of plants in different land biomes (Ren, Chen, and An 2017). We found that the date of the SOS of the forest was the earliest, EOS was the latest, and LOS was the longest; shrubland exhibited the opposite (Fig. S5). The onset of phenology was dependent on the current hydrothermal conditions, and the phenology at different phases affected the timing of the phenology in other periods (Fig. S6). Therefore, this intricate plant-climate regulation mechanism may also indirectly affect the response of land surface phenology to drought stress (Kang, Wang, and Liu 2018;Wang et al. 2020). Different vegetations develop diverse survival strategies and plant-environment feedback systems in the long-term succession process to accommodate frequent drought stress (Takahashi et al. 2020;Gupta, Rico-Medina, and Cano-Delgado 2020). Consequently, land surface phenology exhibited different response intensities and cumulative months when drought occurs. The community of woody plants with deep roots is often highly dependent on mycorrhiza, which allows them to contend for plentiful resources and succeed in the fierce biotope competition . Moreover, the strong deep soil water uptake capacity makes them less vulnerable to surface water supply (Wu et al. 2018;Fan et al. 2017). Conversely, the shallow root system of herbaceous plants can effectively explore soil nutrients, reduce the reliance on symbiotic mycorrhizal fungi and enable rapid responses to periodic resource scarcity and adverse environmental stresses (Xi et al. 2018;Dorji et al. 2013). Furthermore, the water storage capacity and transportation process of plant tissues for woody plants are more complicated than those of herbaceous plants and may help enhance their drought tolerance (Tian et al. 2018). Drought has a profound effect on the growth of trees that occurs during the period of cambium development, and water deficiency can prompt early senescence of leaves, especially old leaves, to minimize the water deficit and redistribute nutrients (Piao et al. 2019). Moreover, the gene-level changes caused by drought stress may induce strong drought resistance to prevent cell death and result in a shift of photosynthetic sources/sinks between different tissues and organs of the plant, thus delaying and extending the EOS and LOS (Xie, Wang, and Silander 2015). We found that the SOS and EOS of the forest ecosystem had the shortest and longest time scales in response to the cumulative effect of drought, respectively. The monospecific structure of arable ecosystems makes them less stable in terms of resistance and resilience ). However, manual intervention and other tillage measures (e.g. autumn harvest and winter irrigation) can enable the onset of land surface phenology to quickly adapt to low water balance.
There was an essential relationship between water availability and land surface phenology of terrestrial ecosystems in this study, the correlation between the phenology and cumulative SPEI varied considerably across water gradients ( Figure 5). The proportion of significant correlations of the cumulative effect on all land surface phenology was the lowest in HU, suggesting that the cumulative effect has the least impact on the surface phenology of the humid area. Contrary to the cumulative effect on SOS, most EOS and LOS were dominantly positively correlated. The trend of r max-cml in the different water gradients was consistent, showing a gradual decrease in the integral correlation from AR to HU, indicating that the cumulative effect on land surface phenology in AR and SAR was higher than in SHU and HU, and suggesting the importance of the cumulative water balance in alleviating the impact of drought on land biomes. The corresponding cumulative months also varied greatly along the water gradients. The cumulative months of SOS response to drought showed a trend of increasing and then decreasing from AR to HU, with the highest in SHU (6.55 months) and the lowest in AR (6.37 months) (Figure 5b). For EOS, SAR was affected by the shortest time scale (6.80 months), and the other regions showed an insignificant difference, with SHU (4.93 months) being the lowest. The cumulative months of EOS showed a downward and then upward trend, the shortest was SAR, and the longest was HU (Figure 5d). The time scales for LOS showed an apparent upward trend from AR to HU, with the shortest in AR (5.77 months) and the longest in HU (7.00 months) (Figure 5f). Generally, the impact intensity of surface phenology in arid and semi-arid areas was higher than that in humid and semi-humid areas, and the cumulative months in humid and semi-humid areas were longer than those in arid and semi-arid areas. This result suggests that the s phenology in SHU and HU is less regulated by the cumulative water balance. Water scarcity at the onset of the phenology period exerts a strong effect on land surface phenology in arid regions (i.e. AR and SAR), as a result of limited water availability, plant growth and development are severely limited in these areas (Wu et al. 2015). Conversely, biomes with less severe abiotic stand conditions (weaker drought stress, i.e. HU and SHU) are less affected by cumulative water availability when land surface phenology occurs. This finding is plausible because plant communities in arid areas are extremely water-sensitive to the amount of available accumulated water, and high-intensity droughts in these areas usually exacerbate water scarcity, thereby limiting vegetation activity (Dorman et al. 2013;Araya et al. 2016), and resulting in a strong drought-phenology correlations. Thus, the r max-cml of phenology-SPEI decreases significantly from AR to HU. In addition, land surface phenology in arid regions reacts to SPEI on a short-term scale, perhaps owing to their rapid response to the restricted moisture available in water-scarce areas; thus, drought-induced damage is mitigated through a wide range of physical acclimation and functional

Responses of Land Surface Phenology to Drought at Lagged Time Scales
The influence of time-lag effect on land surface phenology was determined based on the distribution of r max-lag and the positive/negative correlation ratios between different surface phenology parameters and lagged SPEI (1-month) at different lagged time scales ( Figure 6). Overall, the negative correlation between SOS and 1-month SPEI (53.28%) was slightly larger than the positive correlation (46.72%) (Figure 6a). The distribution of the time-lag effect in different lagged months was concentrated in the short lagged time scales (1-5 months), while the largest was observed in 4 months with a positive in 4.23% and a negative in 4.74% (Figure 6b). The largest positive r max-lag was observed in 1-month (r = 0.428), and the absolute largest negative r max-lag was at 5-month (r = −0.434). The correlation between lagged SPEI and EOS varied considerably across lagged time scales (Figure 6c). The lagged months corresponding to the r max-lag occurrence between EOS and lagged SPEI were mainly distributed 1-3 months, accompanied by a maximum occurring in 1-month (positive in 5.66% and negative in 4.51%). The largest positive r max-lag was observed at 1-month (r = 0.440), while the absolute largest negative r max-lag was at 4-month (r = −0.431) (Figure 6d). The proportion of EOS showing a positive correlation (55.13%) with lagged SPEI was higher than the proportion of negative correlation (44.87%) on all lagged time scales. The drought time-lag effect for LOS was investigated, and the correlation was mainly concentrated in long-term scales of 9-11 months (Figure 6e). The largest proportion was observed in 10 months with a positive of 6.16% and a negative of 2.98%. Generally, the positive correlation between LOS and 1-month SPEI (57.74%) was larger than the negative correlation (42.26%). The largest positive r max-lag was observed in 8-month (r = 0.444). For the negative r max-lag , the largest absolute value was 0.434 in 1-month (Figure 6f).

Spatial Distributions of Drought Time-lag Effect on Land Surface Phenology
The Spatial allocation of r max-lag between different surface phenology parameters and lagged SPEI (1-month) and their corresponding time scales shows the impact of time-lag effects on surface phenology at the global scale (Figure 7). Generally, approximately 46.12% of the vegetated area showed r max-lag , and its distribution was spatially heterogeneous (Figure 7a). The positive-correlation regions (21.55%) were mainly distributed in northern Russia, southern Asia, northern Australia, a small part of western Africa, southwestern North America, and northwestern South America. The negatively correlated regions (24.57%) were concentrated in central Eurasia, southeastern Australia, southeastern Africa, western and northeastern North America, and central South America. In consideration of the lagged time scales, approximately 43.27% of the significant correlation areas were affected by the short-term lagged months (i.e. 1-5 months) (Figure 7b). The spatial patterns with r max-lag and the corresponding lagged months for each pixel between EOS and 1-month SPEI was concluded. Overall, 47.93% of the vegetated area showed r max-lag between the EOS and the lagged SPEI ( Figure 7c). The positively correlated areas (26.42%) were mainly concentrated in most areas of northern Russia, southwest of Eurasia, the eastern part of Australia, and western and northeast parts of South America. The negatively correlated areas (21.51%) were scattered in a small part of northern Russia, central Asia, western Australia, central Africa, northern and southwestern North America, and southeastern South America. A large variation was observed in the spatial distribution of the EOS response's lagged months. The peaks of the lagged time scales were mainly concentrated in 1-3 months (28.67%) (Figure 7d). The LOS of the 50.45% vegetated area was significantly affected by the time-lag effect (Figure 7e). The regions with significant positive and negative correlations accounted for 29.13% and 21.32%, respectively. The positive-correlation regions were mainly distributed in the northeast of Europe, most of northern Russia, eastern Australia, northern and southwestern Africa, eastern and northwestern North America, and the central part of South America. The negatively correlated regions were mainly located in the central Eurasian continent, southwestern Europe, a small part of northern Russia, northeastern China, Western Australia, and most of northern North America. The peak lagged months corresponding to the occurrence of r max-lag were mainly concentrated on the long-term scale, that is, 9-11 months, accounting for 26.83% of the significant-correlation area (Figure 7f).

Drought Time-lag Effect on Land Surface Phenology at Different Land Biomes and Water Gradients
Biological implication suggests that parallel stages of vegetation life history are influenced by water availability, and that water accumulation values in the previous phase of the plant growth and development indirectly influence the next phase of vegetation growth and development. Therefore, water accumulation in the preceding stages of vegetation growth affects land surface phenology through a time-lag effect (Cong, Shen, and Piao 2017). The time-lag effect on land surface phenology varied substantially among the different land biomes (Figure 8). SOS was mainly negatively correlated with 1-month SPEI, whereas EOS and LOS were positive correlation with lagged SPEI, and the positive correlation took up a large area. The distribution of r max-lag showed a common trend in the different land biomes, that is, it increased from forest to cropland and then decreased. The largest r max-lag was observed in shrubland, followed by grassland, forest, and cropland. The corresponding lagged months of the different phenology varied greatly. The lagged months of SOS in grassland were the shortest (6.35 months), and the difference was not significant for the other vegetation types. The EOS in the forest had the longest lagged months (6.61 months), followed by 6.42 months in cropland and 6.24 months in shrubland and grassland. For LOS, the lagged months showed an upward trend from forest to cropland, with the longest being cropland (6.71 months) and the shortest being forest (6.55 months). Overall, the LOS of the different land biomes had the largest proportion (47.62%) affected by the time-lag effect, followed by EOS (45.08%) and SOS (43.36%). LOS had the longest response time scale of the time-lag effect (6.61 months), followed by SOS (6.45 months) and EOS (6.38 months).
The time-lag effect on land surface phenology across water gradients was also explored. The r max-lag and corresponding lagged months varied greatly across different water gradients (Figure 9). The area where SOS was significantly affected by the negative correlation of the lagged month (53.38%) was slightly larger than the area affected by the positive correlation (48.11%). EOS and LOS showed the opposite result, which was mainly positively correlated with 1-month SPEI. The r max-lag showed an apparent downward trend from AR to HU, where AR was the largest and HU was the smallest. The lagged months of SOS were the longest in HU (6.61 months), followed by AR (6.31 months), SHU (6.26 months), and SAR (6.16 months). For EOS, the lagged months did not vary significantly in the different water gradients, with the longest being SAR (6.34 months) and the shortest being SHU (6.31 months). The lagged time scale of LOS in HU (6.77 months) was greater than in other regions. Overall, LOS in the different water gradients had the largest proportion (46.58%), followed by EOS (44.26%) and SOS (42.51%). The longest time scale in response to the drought time-lag effect was observed for LOS (6.51 months), followed by SOS (6.33 months) and EOS (6.31 months). Generally, there were similar patterns in drought time-lag and cumulative effects on land surface phenology at both land biome and water gradient levels. However, the lagged time scales were longer than those of the cumulative time scales. The analysis of the above results confirms our hypothesis that terrestrial biological communities are influenced by the availability of water at the onset of phenology and are also accompanied by a time-lag effect of cumulative available water in the months before phenology.

3 Drivers of Land Surface Phenology in Response to Drought Stress
Many studies have established that the increase in evapotranspiration exacerbates the intensity of drought not only in areas of precipitation scarcity but also induces drought in areas with abundant rainfall (Cook et al. 2014;Zhao et al. 2021). Even though SPEI only depends theoretically on precipitation and evapotranspiration, how these two variables contribute to drought is complex and varied (Vicente-Serrano et al. 2015). Consequently, different hydrothermal factors affect the occurrence time and frequency of drought and thus indirectly influence the SPEI-phenology relationship. Therefore, we calculated the timing of phenology in different land biomes along the water gradient, combined it with hydrothermal conditions and plant physiological factors, and explored the regulation of the complicated vegetation-environment relationship in the effect of drought on phenology. The Pearson correlation heatmap ( Figure 10).
And the RDA plot ( Figure 11). Showed the potential impact of hydrothermal conditions and plant physiological factors on the cumulative and time-lag effects of land surface phenology in response to drought. The RDA plot showed that the first and second axes accounted for 47.77% and 24.14%, respectively, of the variance of the cumulative effect on SOS and its influencing factors (Figure 11a). When combined with Pearson correlation analysis, the r max-cml between SOS and cumulative SPEI (SCR) was a significantly positive correlation (p < 0.05) with evaporation (e). The corresponding cumulative months (SCS) showed an extremely significant positive correlation with ROG (p < 0.01), and significantly positive correlation with SOS and ROS, an extremely significantly negatively correlated with t2m, and a significant negative correlation with ssr. Regarding the cumulative effect on EOS, the explanations of the first and second axes of the RDA diagram accounted for 42.25% and 31.71%, respectively (Figure 11b). The r max-cml between EOS and cumulative SPEI (ECR) showed an extremely significantly positive correlated with e, an extremely significantly negative correlated with LOS, and a significant negative correlation with tp. The corresponding cumulative months (ECS) had an extremely significant positive correlation with str and a significantly positively correlated and significantly negatively correlated with EOS and ssr, respectively. The impact of different forces on the cumulative effect of LOS response to drought had large variations. The first and second axes of the RDA plot accounted for 47.84% and 28.96% of the variance, respectively (Figure 11c). The r max-cml between LOS and cumulative SPEI (LCR) showed an extremely significantly positively correlated with e, an extremely significantly negatively correlated with LOS, and a significant   negative correlation with tp. The corresponding cumulative months (LCS) had an extremely significant positive correlation with ROG and ROS, a significantly positively correlated with str, and an extremely significantly negatively correlated with ssr and t2m. In summary, the drought-phenology relationship was affected not only by environmental conditions but also by the occurrence time and growth rate of phenology. For the time-lag effect, the first and second axes of the RDA chart accounted for 40.38% and 33.00%, respectively, of the variance of the time-lag effect on SOS and its driving factors (Figure 11d). The r max-lag between SOS and 1-month SPEI (SLR) was extremely significantly positively correlated with e (p < 0.01), extremely significantly negative correlation with EOS and tp, and significantly negative correlation with LOS. The corresponding lagged months presented a significantly positively correlated with str (p < 0.05). The first and second axes of the RDA chart could explain 42.96% and 28.72% of the variance of the relationship between the time-lag effect of EOS and its influencing factors, respectively (Figure 11e). The r max-lag between EOS and 1-month SPEI (ELR) had an extremely significant positive correlation with e. The corresponding lagged months (ELS) presented a significantly positively correlated with tp and a significantly negatively correlated with SOS. In the RDA diagram of the impact of different factors on the LOS response to the drought time-lag effect, RDA1 and RDA2 had values of 39.36% and 26.80%, respectively ( Figure 11f). r max-lag (LLR) and e showed a significantly positive correlation, but the time scale (LLS) had no significant correlation.
Overall, the correlation of phenology with SPEI (i.e. the degree of impact by drought) was mainly controlled by hydrological conditions, and the corresponding time scales (i.e. the length of time for phenology to respond to drought events) was mainly related to thermal conditions. Moreover, the onset of phenology and the growth/senescence rate of plants also had a complex impact on the phenology-SPEI relationship. Specifically, evaporation had a positive impact on the correlation between phenology and SPEI, whereas precipitation alleviated the correlation. Net thermal radiation extended the time scales of the phenology response to drought, whereas temperature and net solar radiation did the opposite. The intensity of this regulation in the cumulative effect was higher than that in the time-lag effect, indicating that hydrothermal conditions may shift the occurrence time of phenology by regulating the growth/senescence rate of vegetation and may thus control the phenology-drought interaction. Lastly, this study was concerned only with land surface phenology in response to the drought cumulative and time-lag effects, and the interaction between the two effects was not considered. In addition, human activities (e.g. urban expansion, afforestation) and non-climatic disturbances (e.g. plant community succession, nitrogen, and phosphorus deposition) also impact the vegetation-drought relationship. Future research is possible to combine mathematical models to comprehensively assess how land surface phenology responds to environmental stresses to develop a better understanding of how ecosystem carbon fluxes change, and manage vegetation more effectively, contributing to our understanding of vegetation-environment interactions.
Overall, the correlation of phenology with SPEI (i.e. the degree of impact by drought) was mainly controlled by hydrological conditions, and the corresponding time scales (i.e. the length of time for phenology to respond to drought events) was mainly related to thermal conditions. Moreover, the onset of phenology and the growth/senescence rate of plants also had a complex impact on the phenology-SPEI relationship. Specifically, evaporation had a positive impact on the correlation between phenology and SPEI, whereas precipitation alleviated the correlation. Net thermal radiation extended the time scales of the phenology response to drought, whereas temperature and net solar radiation did the opposite. The intensity of this regulation in the cumulative effect was higher than that in the time-lag effect, indicating that hydrothermal conditions may shift the occurrence time of phenology by regulating the growth/ senescence rate of vegetation and may thus control the phenology-drought interaction. Lastly, this study was concerned only with land surface phenology in response to the drought cumulative and time-lag effects, and the interaction between the two effects was not considered. In addition, human activities (e.g. urban expansion, afforestation) and non-climatic disturbances (e.g. plant community succession, nitrogen, and phosphorus deposition) also impact the vegetation-drought relationship (Aronson, Goulden, and Allison 2019;Ferlan et al. 2016). Future research is possible to combine mathematical models to comprehensively assess how land surface phenology responds to environmental stresses to develop a better understanding of how ecosystem carbon fluxes change, and manage vegetation more effectively, contributing to our understanding of vegetation-environment interactions.

Conclusion
The drought cumulative and time-lag effects for global land surface phenology from 1982 to 2015 were identified according to the maximum correlation between different phenological parameters and SPEI (1-to 12-month and 1-month). The principal conclusions were summarized as follows: (1) The cumulative effect on SOS, EOS, and LOS accounted for 25.36%, 26.43%, and 26.57% of the vegetated area, respectively. Meanwhile, the time-lag effect had a significant effect on SOS in 46.12% of vegetated areas, 47.93% of EOS, and 50.45% of LOS. (2) The drought cumulative and time-lag effects were similar across the different land biomes and water gradients on land surface phenology. Specifically, the r max-cml and r max-lag in shrubland were the largest, the SOS (5.93 months) and LOS (6.55 months) responses to the cumulative and lagged time scales of forest were the shortest, and that of EOS (6.61 months) was the longest. r max-cml and r max-lag were larger in the AR and SAR than in the HU and SHU. The response time scales were longer in the HU and SHU than in the AR and SAR. (3) The correlation between phenology and SPEI was dominated by hydrological conditions, and the time scales were mainly related to thermal conditions. Specifically, evaporation had a positive impact on the phenology-SPEI correlation degree, whereas precipitation alleviated the correlation. Net thermal radiation extended the time scales of the phenology response to drought, whereas temperature and net solar radiation did the opposite. Moreover, the onset of phenology and the growth/senescence rate of plants exerted a complex impact on the phenology-SPEI relationship. This result indicates that hydrothermal conditions may shift the onset of phenology by regulating the growth/senescence rate of vegetation and may thus control the interaction between phenology and drought. Our results indicate that drought stress has significant cumulative and time-lag effects on global land surface phenology. Clarifying these effects would provide valuable information for insights into the variability of land biomes' responses to changing climate and facilitate us to develop vegetation conservation and land degradation reduction strategies more effectively. Future research would combine mechanistic learning models and other accurate predictions of phenology information to further advance our knowledge and understanding of ecosystem carbon fluxes.