Drivers of the pre-season drought thresholds triggering earlier autumn foliar senescence in the Northern Hemisphere

Domain: earth_science
Published: 2025-08-14
Editorial outcome: major_revision

Abstract
Global warming can postpone the autumn date of foliar senescence (DFS). Nevertheless, warming-associated droughts may induce earlier DFS. However, pre-season drought thresholds triggering an earlier DFS (PDT-DFS) are not clearly established. Using site-level DFS data since 1951, satellite-derived DFS data for 1982‒2021, and drought indices, we construct a copula-based Bayesian framework to identify the PDT-DFS over the Northern Hemisphere (>30°N). A higher probability of droughts is associated with an earlier DFS. The DFS for the <10%, <20%, <30%, and <40% quantiles (a lower quantile indicates that DFS will occur earlier under the same drought conditions) results in PDT-DFS values of −2.59, −2.30, −1.80, and −1.63, respectively. The propagation thresholds from meteorological droughts to soil droughts determine the PDT-DFS. However, an increase in resilience and leaf area index hinders the sensitivity of an earlier DFS to droughts. In Sixth Coupled Model Intercomparison Project (CMIP6) simulations, the PDT-DFS increases significantly (p<0.01) under most climatic scenarios in the future. This study provides extensive evidence for the increasing sensitivity of DFS to pre-season droughts and a basis for enhanced predictive responses to such droughts.

Introduction
The autumn date of foliar senescence (DFS, i.e., autumn phenology) plays a crucial role in controlling the length of the growing season in the carbon (C) uptake of terrestrial ecosystems1–4. It is widely acknowledged that global warming delays autumn phenology by reducing the rate of accumulation of chilling units in autumn trees5–8. However, some studies also suggest that global warming has a lower explanatory power for delayed autumn phenology9,10. Warmer autumns may sometimes be associated with an earlier DFS11. This might be caused by the fact that water availability contributes more substantially to interannual variation in DFS than warming12, particularly in arid and semi-arid ecosystems13,14. Under drought conditions, vegetation may modify its strategies for allocating C by prioritizing essential physiological processes over growth, thereby accelerating the senescence of leaves15. Moreover, hydraulic strategies, including adjustments in stomatal conductance and xylem structure, play a crucial role in alleviating water stress; they can also affect the timing of DFS16.

Vegetation is resistant to short-term water deficiencies because of the buffering effect of soil moisture and the fertilization effect of carbon dioxide (CO₂), which can alleviate the impact of droughts on DFS. This resistance safeguards vegetation until it experiences maximum drought stress17. Once this threshold is exceeded, vegetation transitions from highly resistant to highly vulnerable18,19, thus potentially initiating an earlier DFS. However, the spatial distributions and dynamic changes of these critical pre-season drought thresholds triggering an earlier DFS (PDT-DFS) remain unclear despite their significant implications for understanding, simulating, and predicting responses to droughts owing to autumn phenology.

Pre-season droughts related to warming may be crucial causes of an earlier DFS5,6. This process is initiated by the accumulation of abscisic acid and is accompanied by the breakdown of chlorophyll and other pigments, such as β-carotene and lutein, along with the remobilization of nutrients20. This process is also regulated by temperature, solar radiation, water availability, and spring leaf-out dates11,21,22. This suggests that a nonlinear and complex response exists between autumn phenology and drought events.

For example, the compound effects of heat waves and droughts further influence vegetation23 because heat waves can lead to soil and hydrological droughts by propagating meteorological droughts into soil and hydrological systems24,25. This can potentially lead to a more pronounced and earlier DFS, owing to an increase in water stress and the acceleration of physiological responses in vegetation.

However, the asymmetric effects in daytime and nighttime warming may have different impacts on the earlier DFS caused by droughts26. Daytime warming, which generally leads to an increase in the diurnal temperature range, amplifies the adverse effects of drought stress on vegetation growth by decreasing stomatal conductance and photosynthesis27. However, nighttime warming, which generally leads to a lower diurnal temperature range, has a more nuanced relationship with vegetation growth and droughts. For example, studies conducted in growth chambers have observed heightened drought stress during nighttime warming28, whereas field studies have indicated that nighttime warming may reduce the impact of drought stress on vegetation growth29. Thus, the physiological impacts of nighttime warming are still ambiguous. Nighttime warming may elevate respiration rates, which causes a net loss in C over a diurnal period30. However, it might also induce compensatory photosynthesis the following day, thereby enhancing C uptake. It is not clear whether and how pre-season warming and its asymmetric changes intensify changes in the critical drought state (drought thresholds) that result in an earlier DFS.

We construct a framework to identify the PDT-DFS based on copula functions. The advantage of copula functions lies in their ability to capture the nonlinear relationships between pre-season droughts and an earlier DFS, and they can identify drought thresholds that trigger varying degrees of advancement in DFS. We focus on terrestrial ecosystems across the Northern Hemisphere (>30°N) where vegetation dynamics exhibit distinct seasonal variations, with the goal of addressing the following three questions: (1) What are the spatial patterns in the PDT-DFS? (2) How do pre-season heat waves alter the PDT-DFS? and (3) What are the drivers of the PDT-DFS? First, we use copula functions to construct a joint distribution function between pre-season meteorological drought indices (unless otherwise specified, drought thresholds refer to meteorological drought thresholds calculated by the Standard Precipitation-Evapotranspiration Index [SPEI]) and an earlier DFS to identify the joint probability. We then iterate the joint probability distribution to determine the PDT-DFS (see Methods). Secondly, by establishing the joint probability distribution for drought thresholds and heat waves, we evaluate how pre-season heat waves enhance the PDT-DFS (see Methods). We then evaluate the propagation mechanisms among the pre-season heat waves and droughts, and an earlier DFS based on a copula-based Bayesian framework (see Methods). Finally, we analyze the drivers of the PDT-DFS by incorporating the propagation mechanism and other potential drivers into a random forest model (see Methods).

Results
Pre-season drought thresholds triggering an earlier DFS
The results of the copula-based Bayesian framework indicated that an earlier DFS is associated with higher drought thresholds based on both satellite-derived and site-level ground datasets. When the PDT-DFS approaches −0.5, a lower drought intensity is required to trigger an earlier DFS. Conversely, when the PDT-DFS approaches −5.0, a higher drought intensity is needed to trigger an earlier DFS.

For example, the DFS for the <10%, <20%, <30%, and <40% quantiles (a lower quantile indicates that the DFS will occur earlier under the same drought conditions) resulted in PDT-DFS values of −2.59, −2.30, −1.80, and −1.63, respectively (Fig. 1a–d). The deserts and xeric shrubland biomes (DXS) had a higher PDT-DFS compared to other biomes, with PDT-DFS values of −1.75, −1.46, −1.56, and −1.60, respectively, for quantiles <10% to <40% (Fig. 1e). This indicates that the earlier DFS in these biomes is the most sensitive to pre-season droughts.

Moreover, we also found that with the spatial transition of vegetation from isohydric to anisohydric biomes, the PDT-DFS increased (tending toward −0.5) (Supplementary Figs. S1, S2). This indicates that, compared to isohydric biomes characterized by stomatal closure during drought periods, anisohydric biomes can keep their stomata relatively open to maintain photosynthesis, which may exhibit a strong response to droughts.

The DFS based on the site-level ground datasets was less sensitive to droughts compared to the satellite-derived dataset. For example, for the DFS from the quantiles <10%, <20%, <30%, and <40%, the PDT-DFS values were −2.81, −2.58, −2.25, and −1.64, respectively (Fig. 1f). Despite this difference in their sensitivity, the PDT-DFS identified from the satellite-derived data and site-level ground datasets were highly consistent (r=0.92, p<0.01) (Fig. 1g). These results suggest that while satellite-derived data may overestimate the sensitivity of the earlier DFS to droughts, the overall patterns of drought thresholds are consistent between the two datasets.

The PDT-DFS was not identified in regions where droughts are unlikely to be the dominant factor controlling earlier DFS (joint probability between drought and earlier DFS<0.5), or where more intense drought conditions (i.e., SPEI<−5) may be necessary to trigger an earlier DFS. For the DFS from quantiles <10% to <40%, the PDT-DFS were identified in 55.15%, 58.51%, 64.28%, and 66.20% of the regions, respectively (Fig. 1a–d). In central North America, Mediterranean regions of Europe, and north-eastern regions of China, the PDT-DFS was higher (closer to −0.5), which indicated that the earlier DFS in these regions is the most sensitive to pre-season droughts.

In contrast, the PDT-DFS were not observed in most regions of northwestern North America, the central northern part of Siberia, and northern East Siberia (Fig. 1a–d). Similarly, heat wave thresholds triggering an earlier DFS were also not observed in most of these regions (Supplementary Figs. S6, S7a–d). However, in areas where the drought thresholds caused an earlier DFS, the daytime heat wave thresholds were a contributing factor in 66.8–84.7% of the regions (Supplementary Fig. S8), and the nighttime heat wave thresholds were also influential in 61.1–78.4% of the regions (Supplementary Fig. S9). This suggests a potential association between pre-season drought thresholds and pre-season heat wave thresholds in triggering an earlier DFS.

Pre-season heat waves increase the drought thresholds triggering an earlier DFS
A correlation analysis indicated that pre-season heat waves amplified the PDT-DFS based on both satellite-derived and site-level ground datasets (Fig. 2c and Supplementary Fig. S10). However, daytime heat waves had a stronger amplifying effect on the PDT-DFS compared to nighttime heat waves. For example, for DFS from the <10% to <40% quantiles, the correlation coefficients between the drought thresholds and the compound event probability calculated by the standardized maximum temperature index (STImax) and SPEI were 0.76 (p<0.01), 0.74 (p<0.01), 0.70 (p<0.01), and 0.62 (p<0.01) (Fig. 2a, c). In contrast, such correlation coefficients between the drought thresholds and the compound event probability calculated by the standardized minimum temperature index (STImin) and SPEI were 0.64 (p<0.01), 0.62 (p<0.01), 0.59 (p<0.01), and 0.51 (p<0.01) (Fig. 2b, c).

Mechanisms of pre-season drought thresholds triggering an earlier DFS
To explore how the pre-season heat waves amplify the PDT-DFS, we established propagation mechanisms among the pre-season heat waves, pre-season droughts, and earlier DFS based on a copula-based Bayesian framework (see Methods). Pre-season heat waves mediate droughts and thereby lead to an earlier DFS, which is primarily achieved through 14 potential propagation pathways (Fig. 3e). They are divided into three categories. First, pre-season heat waves trigger an earlier DFS by propagating into meteorological droughts; second, pre-season heat waves directly induce soil and hydrological droughts, which cause an earlier DFS; and third, pre-season heat waves propagate into meteorological droughts, which then trigger soil or hydrological droughts, thus resulting in an earlier DFS.

Comparing the propagation thresholds of 14 potential propagation pathways showed that an earlier DFS was the most sensitive to the pathway in which the pre-season heat waves cause meteorological droughts, thus subsequently causing soil droughts. Satellite-derived and site-level ground datasets both showed that thresholds for this propagation pathway were closer to −0.5 (Fig. 3a–d and Supplementary Figs. S11–S21).

Furthermore, we found that the PDT-DFS, identified using the standardized runoff index (SRI) and total water storage anomaly-drought severity index (TWSA-DSI) based on both satellite-derived and site-level ground datasets, were more dependent on daytime heat waves (Fig. 2c and Supplementary Fig. S8). Conversely, the PDT-DFS, identified using the standardized soil index (SSI) derived from both satellite-derived and site-level ground datasets, depended less on daytime heat waves (Fig. 2c and Supplementary Fig. S10). The results of the propagation mechanisms also showed that thresholds for soil droughts (measured by SSI) triggered by daytime heat waves were generally higher than those triggered by nighttime heat waves (Fig. 3a–d and Supplementary Figs. S22–S24). This indicates that respiration may adapt to a higher STImin and partially offset C losses caused by nighttime heat waves31,32.

Drivers of the pre-season drought thresholds triggering an earlier DFS
To further analyze PDT-DFS drivers, we incorporated 72 potential drivers, including 14 propagation mechanism drivers and 58 environmental drivers (Supplementary Table S3), into the random forest model (see Methods). Recursive feature elimination methods were used to select a subset of 16 influential drivers, which explained 0.81–0.86 of the spatial variation in the PDT-DFS (see Methods). The drivers primarily included five categories: propagation mechanism drivers, climatic drivers, vegetation drivers, hydraulic characteristics of vegetation, and topographic drivers.

The PDT-DFS were predominantly governed by propagation thresholds from the SPEI to SSI, which accounted for more than 0.55 of their importance (average values of 100 cross-validations) (Fig. 4a–d). Secondly, the PDT-DFS were governed by propagation thresholds from the SPEI to TWSA-DSI, which accounted for more than 0.19 of their importance (Fig. 4a–d). Partial dependence plots further revealed that the PDT-DFS also increased as propagation thresholds from the SPEI to SSI and TWSA-DSI increased (close to −0.5), (Fig. 4e, f). This suggests that stronger mediation of soil and hydrological droughts by meteorological droughts results in higher drought thresholds (close to -0.5) triggering an earlier DFS.

Furthermore, after excluding the 14 propagation mechanism drivers in the random forest model, the recursive feature elimination method was applied to select a subset of 12 influential predictor drivers from 58 environmental drivers (Supplementary Table S3). These 12 drivers explained 0.72–0.79 of the spatial variation in the PDT-DFS, which were predominantly governed by resilience, and accounted for more than 0.69 of the importance (Supplementary Fig. S25a–d and Supplementary Fig. S26a–d). As resilience increased, the Shapley values (SHAP) tended toward lower negative values, which indicated that a higher ecosystem resilience would hinder the sensitivity of DFS to droughts (Supplementary Figs. S25e and S26e).

Resilience and leaf area index (LAI) served as the primary mechanisms buffering ecosystem sensitivity to droughts based on structural equation modeling (SEM) results. This finding aligned with the results derived from the random forest model after removing the 14 propagation mechanism drivers. For example, we found positive effects of precipitation (PPT, path coefficients ranging from 0.54 to 0.55) and aridity index (AI; higher AI values indicate more humid conditions with path coefficients ranging from 0.77 to 0.78) on resilience that subsequently delayed the occurrence of an earlier DFS (path coefficient approximately −0.05 to −0.06) (Supplementary Fig. S28). Increased precipitation positively enhanced LAI (path coefficients 0.54–0.56), which in turn significantly delayed an earlier DFS (path coefficient around −0.02 to −0.03). These findings indicate that sufficient precipitation enhances vegetation growth potential and effectively mitigates ecosystem sensitivities to droughts (Supplementary Fig. S28).

Some variation existed in the buffering effects of resilience and LAI on ecosystem sensitivity to droughts across different biomes. For example, in temperate broadleaf mixed forests (TBMF), AI notably delayed earlier DFS through enhanced LAI (Supplementary Figs. S29–S30). In contrast, Boreal forests/taiga (BF) and tundra ecosystems (TUN) showed stronger delayed earlier DFS through enhanced resilience. Moreover, in deserts and xeric shrublands (DXS), both AI and PPT significantly delayed earlier DFS through enhanced resilience and LAI, thus reflecting the crucial role of water availability in drought-prone ecosystems.

Projected future changes in pre-season drought thresholds triggering an earlier DFS
Based on the Inter-Sectoral Impact Model Intercomparison Project Phase 3b (ISIMIP3b) datasets from the averages of five global climate models (GCMs) (for the dataset details, see Supplementary Notes S1 and Supplementary Tables S5–S6), we found that the PDT-DFS increased significantly (p<0.01) (Fig. 5a–d) under most of the climatic scenarios over the Northern Hemisphere. This indicates that there will be an increased sensitivity of DFS to pre-season droughts in the future. Compared with the estimates for 1982–2021, more than 51.09% of the areas had higher drought thresholds in 2061–2100 (Supplementary Fig. S31). This indicates that the future DFS would be more sensitive to pre-season droughts in most areas of the Northern Hemisphere. Higher drought thresholds were primarily concentrated in the northern part of the Eurasian continent under the three climate scenarios of SSP1-2.6, SSP3-7.0, and SSP5-8.5 (Supplementary Fig. S31).

Moreover, correlation coefficients between the compound drought and heat wave events (CDHWs) and the PDT-DFS tended to increase under most of the climatic scenarios (Supplementary Notes S3 and Supplementary Figs. S38–S39). Furthermore, correlation coefficients between the sensitivity of vegetation to extreme climate events, including heat waves, droughts, and CDHWs, and the PDT-DFS also tended to increase under most of the climatic scenarios (Supplementary Notes S3 and Supplementary Figs. S40–S44). These results highlight the roles of CDHWs in enhancing PDT-DFS trends. However, the adaptability of vegetation to extreme climates, including heat waves, droughts, and CDHWs, hindered increases in PDT-DFS.

Methods
Satellite-derived DFS datasets
The satellite-derived DFS datasets for 1982–2021 were developed using the normalized difference vegetation index (NDVI). The NDVI data were sourced from the third-generation V1.2 (GIMMS-3G+) dataset, with a 0.0833°× 0.0833° spatial resolution at the global level from 1982 to 2022. NDVI values are reported twice per month. This dataset was assembled from various AVHRR sensors and accounted for various problematic issues, such as calibration loss, orbital drift, and volcanic eruptions.

First, we excluded areas with a multi-year average NDVI of <0.1.
Secondly, we excluded farmland and areas that lacked vegetation based on the moderate-resolution imaging spectrometer Land Cover Version 6 data product (MCD12C1 v006). Third, we applied the Savitzky-Golay smoothing method to the NDVI time series data to reduce noise, while preserving seasonal trends in vegetation dynamics55. In particular, we used a 5×5 moving window, where the window size represented a temporal span of five observations that were centered on each point, with a polynomial order of 3 for the optimal fit. This approach has been widely used in phenological studies because it effectively reduces high-frequency noise, while retaining the signal of vegetation growth and patterns of senescence34,55.

Based on the reconstructed NDVI time series, we extracted the DFS dataset for 1982–2021 over the Northern Hemisphere (>30°N) by applying the Hants-Maximum56, Polyfit-Maximum57, double logistic58, and piecewise logistic methods59. These methods utilize different mathematical models to describe the NDVI time series, which enabled us to identify the critical transition point where the vegetation activity shifts from peak growth into the rapid decline phase, which denotes the acceleration of senescence (Supplementary Notes S2). To reduce the uncertainty of DFS extracted from a single method, we calculated the annual arithmetic average DFS based on these four methods34,60 (Supplementary Fig. S47a).

Additionally, we chose the Mid Greendown band from MCD12Q2, which represents the date when EVI2 last crossed 50% of the segment EVI2 amplitude (MODIS-derived DFS). Since the MODIS-derived DFS reflects the acceleration phase of senescence, it serves as a valuable indicator to validate the consistency and robustness of our DFS dataset. We found that the DFS derived from AVHRR GIMMS-3G+ and MODIS-derived DFS are highly correlated (r=0.73, p<0.01) (Supplementary Fig. S49). Furthermore, we found that calculating the PDT-DFS using AVHRR GIMMS-3G+-derived DFS and MODIS-derived DFS results in similar values (Supplementary Figs. S50, S51). This indicated that the PDT-DFS driven by AVHRR GIMMS-3G+-derived DFS and drought indices were highly robust.

Site-level ground DFS datasets
Site-level ground DFS datasets were sourced from the Pan European Phenology project (PEP725). The site-level ground DFS datasets for PEP725 were based on the autumnal foliar senescence (50%) (BBCH= 94). We selected four temperate tree species, including Aesculus hippocastanum, Fagus sylvatica, Betula pendula, and Quercus robur, and each site had at least 30 years of continuous observational data. Moreover, based on the median absolute deviation (MAD) method, we excluded outliers that deviated from the median by more than 2.5 times61. The results included 283, 203 sites of DFS recorded from 1951–2015 (Supplementary Fig. S52c). The mean number of sites for 1951–2015 for each species was 1203, 1156, 1001, and 996 (Supplementary Table S7).

Drought and heat wave indicators
We used the SPEI function in R to calculate the monthly SPEI values at 1- to 12-month scales using the potential evapotranspiration (PET, mm) and precipitation (Pre, mm) for 1951–2021. The monthly PET and Pre were sourced from CRU-TS4.08 of the Climatic Research Unit Time Series (CRU) with a spatial resolution of 0.5°.

The monthly soil moisture values (m3·m−3) for 1951–2021 were sourced from the ERA5-Land reanalysis dataset. The soil moisture data consisted of four layers (0–7cm, 7–28cm, 28–100cm, and 100–289cm) with an original spatial resolution of 0.1°×0.1°. We used a weighted average method to calculate the soil moisture for the top 100cm in depth based on the first three layers62. Based on the methods for equiprobability transformation, we standardized the soil moisture as SSI at 1- to 12-month scales63.

Terrestrial water storage anomalies (TWSA) were obtained from Yin et al.64. These data provide a longer time series compared with the terrestrial water storage measured by the Gravity Recovery and Climate Experiment (GRACE) satellite mission. TWSA-DSI was defined as follows:

TWSA-DSI = (TWSA_i,j - TWSA_j) / δ_j (1)

where SPEI and SSI indicate meteorological droughts and soil droughts, respectively. We acknowledge that soil droughts are not a formal drought type. However, they are often used to represent agricultural droughts. We referred to the SSI as soil droughts rather than agricultural droughts because our analysis excluded areas affected by humans, such as farmlands and urban zones. We also used SRI and TWSA-DSI to represent hydrological droughts (Supplementary Table S8).

Furthermore, we collected monthly maximum temperature (Tmax, °C) and monthly minimum temperature (Tmin, °C) data from CRU-TS4.08 with a spatial resolution of 0.5°. Based on the methods for equiprobability transformation, we first converted Tmax and Tmin to standardized temperature indices at 1–12-month scales (STImax, and STImin, respectively). We then multiplied STImax, and STImin by -1 to conform to the same direction as the drought indices; i.e., larger negative values indicate more severe heat waves (Supplementary Table S8).

To match the spatial resolution of the DFS, we used an area-weighted mean approach to adjust the drought indices to a 0.083°× 0.083° spatial resolution65. Moreover, to verify the impact of the differences in resolution on the drought thresholds that trigger a DFS, we also resampled the climate indices to a 0.5°× 0.5° spatial resolution using an area-weighted mean approach. We generated spatial distribution maps of the PDT-DFS at 0.083° × 0.083° and 0.5° × 0.5°, respectively (Supplementary Figs. S53–S56 and Tables S9–S12). We found that the PDT-DFS is not sensitive to resolution (Supplementary Notes S5).

Auxiliary dataset
The auxiliary dataset included global land use datasets (Supplementary Fig. S52b and Table S2) and global biomes (Supplementary Fig. S52a and Table S1). Based on the MCD12C1 v006 dataset, land use pixels that remained constant between 2001 and 2020 were analyzed, while excluding regions classified as barren, terrain, ice formations, and urban areas. The Earth was divided into eight biomes based on the Terrestrial Ecoregions of the World (TEOW) dataset66. We analyzed the PDT-DFS in different biomes.

Quantifying the joint probability distribution of the pre-season drought indices and an earlier DFS
First, we used copula functions to construct a joint distribution function between the pre-season drought indices and DFS. The copula functions can overcome the shortcomings of counting the co-occurrence rate between pre-season drought indices and an earlier DFS with few samples44,67.

We performed linear detrending of the DFS dataset. The detrended DFS dataset was then converted to a percentile data set to represent varying degrees of the advancement of DFS68. Lower percentiles of DFS represent earlier advancements of DFS, while higher percentiles of DFS indicate delayed DFS. We used DFS≤10th percentile, DFS≤20th percentile, DFS≤30th percentile, and DFS≤40th percentile to represent the advancement of DFS.

To account for the pre-season drought legacy effects on DFS34,69, we identified the optimal response months of DFS to pre-season drought indices (SPEI, SSI, SRI, and TWSA-DSI) by calculating the absolute correlation coefficients between the DFS and drought indices for each pre-season month (1–12) during the period from 1982 to 2021 (Supplementary Figs. S57, S58 (a)–(d)).

Furthermore, we aggregated daily DFS values into months. The optimal response months of DFS to pre-season drought indices for each pixel were then calculated based on the 1–12 months that preceded the DFS month (Supplementary Fig. S47 (d), (e)). For example, if a DFS pixel is located in August, the length of its pre-season 1–12 period will start from July and move backwards.

Moreover, owing to the limited efficiency of parametric distributions caused by uncertainties in the estimation of parameters and their inapplicability to random variables <070, a kernel distribution was utilized to fit the marginal distributions and in subsequent calculations24. A copula function was used to construct the joint probability distribution between the DFS (x) and SPEI (y) as Formula (2) where F_X and F_Y are the marginal distributions of x, y. C is the cumulative distribution function of copula as shown below:

F_XY(X,Y) = C(F_X(X), F_Y(Y)) (2)

The kernel density estimator is defined by Formula (3) where x, random samples; n, sample size; K, kernel smoothing function; h, bandwidth as shown below68.

^f(x) = 1/(nh) * Σ_{i=1}^n K((x-x_i)/h) (3)

The bivariate copula function was applied to establish the joint distribution function of DFS and pre-season drought indices. Taking SPEI as an example, it was described as Formula (4) shown below:

F_{DFS|SPEI}(DFS≤dfs|spei_u1 ≤SPEI≤spei_u2) = F(DFS≤dfs, spei_u1 ≤SPEI≤spei_u2) / F(spei_u1 ≤SPEI≤spei_u2) = (F_{SPEI,DFS}(spei_u2, dfs) - F_{SPEI,DFS}(spei_u1, dfs)) / (F_{SPEI}(spei_u2) - F_{SPEI}(spei_u1)) (4)

where F_{DFS|SPEI}(DFS≤dfs|SPEI≤spei) represents the joint cumulative probability distribution function (CDF) between DFS below a certain threshold (DFS≤dfs) and SPEI within a certain range (spei_u1 ≤SPEI≤spei_u2). The ranges of SPEI are defined as follows when calculating the joint distribution (all drought, [−0.5, −5), mild drought, [−0.5, −1), moderate drought [−1, −1.5), severe drought [−1.5, −2), extreme drought [−2, −5)). F represents the CDF.

The specific implementation of the bivariate copula function for DFS and pre-season drought indices was as follows
First, we performed normal fitting for the DFS and pre-season drought indices. For the DFS and, SPEI, SSI, and TWSA-DSI, we used the normal distribution and generalized extreme value distribution for fitting, and then selected the optimal marginal distribution based on the Kolmogorov-Smirnov (K-S) test and root mean square error (RMSE). We then used five copulas—Clayton function, Frank function, Gumbel function, Joe function, and t function—as potential joint distribution functions (F) for the DFS and pre-season drought indices. The optimal fitting function for each grid was selected based on the Akaike Information Content (AIC)25,44. This represents the optimal joint probability distribution between the DFS and different pre-season drought indices.

The joint probability distribution between the DFS and different pre-season drought indices ranges from [0, 1]. A higher joint probability indicates a lower likelihood of the simultaneous occurrence of DFS and droughts, while a lower joint probability indicates a higher likelihood of the simultaneous occurrence of DFS and droughts.

Based on the four scenarios (i.e., earlier DFS—DFS<40th percentile, DFS<30th percentile, DFS<20th percentile, and DFS<10th percentile) and four drought intervals (mild drought, moderate drought, severe drought, and extreme drought), we established the possible spatial distribution characteristics of 16 joint probability distributions between an earlier DFS and droughts. As drought severity increased, the probability of triggering an earlier DFS also increased. This implies that more severe pre-season droughts are more likely to trigger an earlier DFS (Supplementary Figs. S63–S66).

Quantifying the pre-season drought thresholds triggering an earlier DFS
Based on the copula function, we iteratively applied the joint probability distribution between the pre-season drought indices and DFS to establish a copula-based Bayesian framework to identify the PDT-DFS as described in Formula (5) (Fig. 6) shown below:

F_{DFS|SPEI}(DFS≤dfs|spei_{i-0.1} ≤SPEI≤spei_i) = (F_{SPEI,DFS}(spei_i, dfs) - F_{SPEI,DFS}(spei_{i-0.1}, dfs)) / (F_{SPEI}(spei_i) - F_{SPEI}(spei_{i-0.1})) (5)

Based on the optimal response months for the DFS and pre-season drought indices for each grid, we constructed the joint distribution of DFS and pre-season drought indices at the maximum timescale, pixel by pixel. Taking DFS≤40% and pre-season SPEI as an example, the SPEI was iterated from −0.5 with intervals of −0.1 until it equated to −5.0. The joint probability of SPEI and DFS was calculated at each iteration. When the joint probability ≥0.5, the iterative process was terminated, and the corresponding SPEI interval was returned. The left side of the interval was taken as the triggering thresholds, which is the SPEI that corresponds to the induction of earlier DFS. A lower triggering threshold (SPEI close to −5.0) indicates that more extreme drought is needed to trigger an earlier DFS. A higher triggering threshold (SPEI close to −0.5) indicates that a milder drought can trigger an earlier DFS. If the entire iterative process ends, and we do not retrieve the corresponding drought thresholds, then no drought threshold exists that triggers an earlier DFS for that pixel.

For the analysis of each pixel, we utilized a 3×3 moving window to execute a copula-based Bayesian framework to identify the PDT-DFS. In particular, we aggregated the values from the adjacent 3×3 pixels to increase the sample size by a factor of 971. Moreover, we also utilized 1×1 and 5×5 moving windows to execute a copula-based Bayesian framework to identify the PDT-DFS. We found little difference between using 1×1, 3×3, and 5×5 moving windows to execute a copula-based Bayesian framework to identify the PDT-DFS (Supplementary Fig. S67). Additionally, executing the copula-based Bayesian framework with a 3×3 moving window to identify the PDT-DFS resulted in good marginal distribution characteristics, which suggested a robust copula-based Bayesian framework (Supplementary Figs. S68–S75). Considering the trade-off between the computational cost and sample size of executing a copula-based Bayesian framework, we ultimately used a 3×3 moving window to identify the PDT-DFS. Moreover, we also assessed the robustness and uncertainty of the copula-based Bayesian framework (Supplementary Notes S4–S5).

Quantifying the pre-season heat waves that amplify drought thresholds triggering an earlier DFS
We established a joint distribution function for the drought thresholds and corresponding monthly STI using Formula (6) as shown below:

F_{SPEI|STI}(SPEI≤spei|STI≤sti) = F(SPEI≤spei, STI≤sti) / F(STI≤sti) = (F_{SPEI,STI}(spei, sti)) / F_{STI}(sti) (6)

where F_{SPEI|STI}(SPEI≤spei|STI≤sti), CDF of drought (SPEI≤spei) and heat wave (STI≤sti) conditions. SPEI≤spei represents drought thresholds at each location; STI≤sti represents STI ≥ to a certain heat wave threshold. In this study, we selected STI≤-2 to establish the joint probability distribution of pre-season heat waves and drought thresholds. If no PDT-DFS were detected in the region, STI≤-2 and SPEI≤-2 were used to represent the joint probability distribution of pre-season heat waves and drought thresholds for that region. This represented relatively extreme compound heat wave and drought events in the region.

Quantifying the mechanisms of propagation from the pre-season heat waves to an earlier DFS
Previous studies suggest that warming-related droughts, which lead to water resource limitations, may in turn stimulate an earlier DFS9,69. Moreover, meteorological droughts can propagate into soil and hydrological systems24,25,47, thus causing soil and hydrological droughts. Therefore, based on a copula-based Bayesian framework, we established the propagation mechanisms among the pre-season heat waves, pre-season droughts, including SPEI, SSI, SRI, and TWSA-DSI, and an earlier DFS.

Based on Formula (3), we calculated 15 sets of propagation thresholds (Fig. 1a–d, Supplementary Figs. S3–S5a–d, and Supplementary Figs. S11–S21). These included the following: SPEI thresholds triggering an earlier DFS; SSI thresholds triggering an earlier DFS; SRI threshold triggering an earlier DFS; TWSA-DSI thresholds triggering an earlier DFS; STImax thresholds triggering meteorological droughts (SPEI as a proxy for meteorological droughts); STImax thresholds triggering soil droughts (SSI as a proxy for soil droughts); STImax thresholds triggering hydrological droughts (SRI and TWSA-DSI as proxy for hydrological droughts); STImin thresholds triggering meteorological droughts; STImin thresholds triggering soil droughts; STImin thresholds triggering hydrological droughts; SPEI thresholds triggering soil droughts, and SPEI thresholds triggering hydrological droughts.

These form 14 potential pathways that could lead to an earlier DFS, which are divided into three categories (Fig. 3h). We calculated the drought triggering thresholds for an earlier DFS under different scenarios (quantiles <10%, <20%, <30%, and <40%) and different drought indices (SPEI, SSI, SRI, and TWSA-DSI). We also calculated the drought and heat wave triggering thresholds for other types of droughts under different scenarios (quantiles <−0.5, <−1.0, <−1.5, and <−2.0) different drought indices (SPEI, SSI, SRI, and TWSA-DSI). When the average thresholds of a potential pathway are closer to −0.5, it is more likely that the heat waves cause an earlier DFS by influencing droughts along that pathway.

Quantifying the drivers of pre-season drought thresholds triggering an earlier DFS
This study used the random forest method to identify the drivers of the PDT-DFS. A total of 72 potential drivers were selected (Supplementary Table S3), and they are categorized into six main groups, including soil, topography, climate, vegetation, vegetation hydraulic characteristics, and propagation mechanism drivers (14 propagation thresholds). As described by Green et al.72 and Deng et al.73, the recursive feature elimination method was used to select the variable combination with the highest coefficient of determination (R²). The number of drivers used for splitting each node was set to 2, and the model contained 500 trees with the goal of maximizing the R².

SHAP values assessed the contribution of each predictor combination to the outcome. Partial dependence plots of the SHAP values can be used to determine the impact of a single variable on the response, independent of the interference of other drivers. A SHAP value >0 indicates that the driving variable promotes the dependent variable; conversely, a SHAP value <0 indicates that the driving variable obstructs the dependent variable.

Based on the SEM, we estimated and examined causal relationships between PDT-DFS and influential environmental factors. We selected five influential environmental factors to discern their mechanisms influencing the PDT-DFS. These five factors were selected according to the random forest model after removing the 14 propagation mechanism drivers. The relative importance of these five variables was consistently ranked in the top five across different DFS percentiles (<10%, <20%, <30%, and <40%), and their ranking remained unchanged, demonstrating that these indicators have an important and stable influence on the PDT-DFS. We used the lavaan package74 in R 4.2.0 to compute standardized path coefficients for the predefined path diagram, calculating them as the product of the standardized coefficients along each path.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Review of “Preseason heatwaves amplify the drought thresholds triggering earlier autumn foliar senescence over the Northern Hemisphere”.\nThe authors proposed a new approach for quantifying the impacts of heatwaves on drought thresholds, using the Coupla-based Bayesian framework. The authors tried to answer how preseason drought thresholds trigger earlier DFS (sate of foliar senescence), and how heatwaves might amplify the thresholds. Finally, the authors found that there is an increasing trend in the thresholds over more than 56% regions, and they are higher when heatwaves occurred during droughts. The CMIP6 simulations also verified this increase in the thresholds. Generally, I think this is an interesting idea for analysis of drought and heatwave impacts. However, I have some concerns on this work:\n1, the authors showed the difference between satellite and site level results, but not clearly demonstrated what are the causes of them. Site level phenology is different from the NDVI derived phenology both for scales and the definitions. Therefore, the authors need to think about how to properly explain their results and what is the consistence in them.\n2, different drought indices were used and also showed discrepancies, which I would say are as expected, giving they represent different aspects and types of droughts. Therefore, I would ask the authors to explain why you need to choose different indices and what you most want to demonstrate. I am not satisfied with the discussion parts on these differences, as they are not convincing. More mechanistic analysis or explanations are necessary.\n3, the current figures are full of information but without clear clue for what is the key one. I am confused about the too many elements and subplots. I would suggest the authors to refine the results and put the most important ones in the figures.\n4, the impacts of droughts and/or heatwaves on vegetation phenology are the result of physiological responses. I would say that the hydraulic properties are very important and need to be considered in the analysis, while current work is more about the surficial explanations, making the results less convincing to readers.\nSome more detailed comments:\nLine 15: predictions of what?\nLine 21: how the warming contributes to autumn phenology, via which? Needs to be more precise.\nLine 25-26: not carbon allocation strategy? And also hydraulic strategy?\nLine 28: this is quite specific for plant species, and hydraulic properties.\nLine 41: I feel that you need to summarize the co-control of heatwaves and droughts on DFS, as this is one part of your analysis.\nLine 86-87: this is what we expect, right?\nLine 88-89: are you sure these regions are all arid and semi-arid regions?\nLine 91: maybe more temperature control?\nFigure 1: You need to explain what is the inset representing for?\nIf it is the density distribution of spatial data, then I think the x=-2.5 is not correctly marked, as you showed -2.5 is green color in you colorbar...\nLine 136-138: much higher than satellite data, why?\nFigure 2: I feel there is not need to show time series changes in thresholds in the main text. Too much information. You can move them to SI.\nLine 228-229: This is reasonable, as soil water is more linked to vegetation physiology, than meteorological droughts, given root water uptake taking place from soil, especially root zone.\nLine 231-234: This is reasonable, as soil water is more linked to vegetation physiology, than meteorological droughts, given root water uptake taking place from soil, especially root zone.\nLine 242-243: This is not a strong point to me.\nI would expect that the increaseing trend is probably due to the more co-control of temperature and droughts.\nTherefore, it is necessary to conduct analysis on the probability of drought-heatwave in the history and future, in order to answer this question.\nLine 260: ‘study’, since only one reference cited.\nLine 263-264: The asymmetry is also due to the physiological changes during daytime and nighttime. At night, the stomata is nearly closed, which could buffer the droughts and heatwaves impacts.\nI suggest to make more analysis on this mechanisms, instead of citing just one study, which is not convincing.\nLine 300-302: I think the acclimation of plants to droughts and heatwaves is also important in a long time span until the end of century. This is also not considered in the current study.", "reviewer_2": "The manuscript (MS) deals with a highly relevant issue regarding the effect of drought on the autumn date of foliar senescence (DFS) in northern latitude ecosystems. Phenology modelling is a hot topic in Earth System Model development, therefore the presented research has significance, and the results could be interesting for the wider audience. The main concerns about this study are related to the mechanisms underlying drought effects and thresholds, as well as the causality versus correlation between heatwaves and drought thresholds. The analyses and interpretation of future projections need to be careful since the uncertainty of CMIP simulations. Overall, the issues with methods and organization and writhing (including but not limited to language) should be fully addressed before the manuscript can be considered for publication.\n1) One main concern about satellite-based DFS is the mismatch of spatial resolution between NDVI data and climatic data.\nIn the study, the authors resampled the DFS into 0.5 degree to match climatic data. However, this would bring large uncertainty of conclusion, especially the range of drought threshold. Previous studies tended to resample climatic data into the resolution of phenology to avoid this problem. Thus, it is necessary to justify the conclusion that is not sensitive to the mismatch of spatial resolution.\n2) Current phenology studies often tried to predict future trends using climatic projections to make the story more interesting.\nBut the uncertainty of data itself makes the predictions unreliable. In this study, the authors used only one experiment \"r1i1p1f1\" from CMIP6 for future prediction, making results and conclusions unconvincing. For climatic data, ISIMIP data has been corrected, which might be more suitable. Besides, the future DFS was estimated by LAI instead of NDVI.\nLAI represents not only canopy greenness information (similar to NDVI) but also structural characteristics. The usage of LAI in DFS estimate, thus, is not common.\n3) The study identified the drought threshold of DFS from site to global scales, however, did not provide biogeophysical mechanisms underlying the drought effects, which should be the necessary step before investigating drought threshold. In other words, audience should understand why drought could advance or delay DFS before being interested in the drought threshold and its response to heatwaves. The mechanism could be related to plant traits, environmental conditions, and ecosystem resilience. The authors should put more efforts into mechanisms instead of results of statistical analyses.\nReviewed by Chaoyang Wu，IGSNRR", "reviewer_3": "This manuscript identified the drought thresholds triggering earlier DFS over the Northern\nMajor comments\n1. Drought was mainly caused by warming temperature in autumn. How to determine where the DFS variation was driven by drought rather than warming temperature?\n2. DFS was used for autumn date of foliar senescence. However, this no definition was provided regarding the exact means of the satellite-derived autumn date in this study. It should be noted that foliar senesce is phase during vegetation development from the start of senescence to the end of senescence.\n3. DFS was derived from AVHRR GIMMS data at 0.0833° pixel. It is not clear how the NDVI time series was smoothed using the Savitzky-Golay method with a 5 × 5 moving window. Moreover, no information was provided regarding the aggregation from 0.0833° to 0.5°.\n4. The DFS was converted to percentile with DFS<10th, 20th, 30th, and 40th. For AVHRR data from 1982-2021, the number of DFS for a pixel was 4 for DFS<10th and 8 for DFS<20th. How could the joint probability distribution function be generated using such small samples for a pixel? Was the function reliable? It is needed to understand the how good the established join distribution functions.\n5. The variables in the equations were not clearly defined. For example, what was speiu and (speiu1 SPEI ≤ speiu2) in line 401?\n6. The drought index was obtained before the average months (September). This could be reasonable in high latitudes. In the regions from 45° to 30° latitudes, the foliage senesces generally started from early October to November. How about the drought impacts occurred in September (or event early October in southern regions)? Based on Fig. S13, there were various regions with DFS earlier than August. How was the drought index used for processing?\n7. The drought index was calculated using monthly data from 1-12 months. The preseason months used in this study varied with pixels. What was the biophysical reason?\n8. Based on Fig. S13a although the color legend was really bad, it looks like the DFS was early in southern regions and later in northern regions. This pattern is very clear in North America. This reviewer thought the pattern should be opposite. Thus, this review suspected the results from this study.\n9. DFS was a variable in a daily basis. The DFS interannual variation in a pixel should generally less than half months. However, the drought index was generated in a monthly scale. This means the DFS values with a difference less than 15 – 30 days would correspond to the same drought index. How could this mismatch be overcome? or how big the impacts?\n10. To quantify the drought thresholds for earlier DFS, a 15-year sliding window was used. Was it also used the samples with DFS<10th, 20th, 30th, and 40th? If so, how many samples were used to generate the distribution functions? Was the function reliable?\n11. Line 43-46, how about the daytime warming and nighttime warming occur together? This could be the most cases in the real world.\n12. Lin 81-83: “with DFS values of -2.20, -1.94, -1.54, and -1.08, respectively, for quantiles <10% to <40% (see supplemental information Table S5).” However, Table S5 did not show the related information.\n13. Line 86-90: Does this means the central North America and Yangtze River Basin belong to arid and semi-arid ecosystems?\n14. The temporal trend of drought thresholds increased in mid-latitudes while decreased in high latitudes. What was the biophysical explanations?\n15. The discussions mainly summarized the results without deep explanations, particularly the biophysics.", "reviewer_4": "General Comments:\nIn this study, the authors tried to develop a copula-based Bayesian framework for identifying the drought thresholds triggering earlier date of foliar senescence, and the authors also used multiple datasets (site-level data, satellite-derived data, and climate models) to verify the reliability of the results. Copula function is a good method to quantitatively analyze the nonlinear correlation of two or more variables, and Bayesian algorithm can reduce the uncertainty of the framework. Therefore, from the perspective of method application, this idea is feasible. Having read your manuscript, I believe it is within the scope of Nature Communications, but the current version is not ready for publication unless the following comments were answered.\nTherefore, I recommend a major revision for this manuscript, and I hope that the authors benefit from the following specific comments for further improvement.\nSpecific Comments (in line number order)\n1. Line 31-35: The author uses three \"however\" in succession in the introduction, which makes the focus of the purpose of this study unclear and confuses the reader. Please make it clear.\n2. Line 58-59: Between the equator and 30°N, there are many regions rich in vegetation, such as southern China, the mountains of the Indian peninsula, and the east and west coasts of Mexico. Human activities are frequent in these regions, and phenological changes have a significant impact on human society. Why does this study only discuss the region north of 30°N?\n3. Line 90-91: The argumentation of “In contrast, no drought thresholds were detected in the northern regions of central Siberia, indicating stronger drought resistance” is inadequate. No drought thresholds were detected in somewhere does not indicate the stronger drought resistance, and further analysis of the reasons for no drought thresholds were detected is needed.\n4. Line 130-131: Deserts and xeric shrublands biomes showed the most significant increases. Vegetation species in deserts and xeric shrublands are mostly plants with developed roots but degraded leaves (e.g. Caligonum, Haloxylon ammondendron) . Are the date of foliar senescence of these plants representative? Please analyze which region has seen the most significant changes in trends except deserts.\n5. Line 168-172: The results based on TWSA-DSI and SSI are opposite. Were the reasons further analyzed in this study?\n6. Line 189-201: Please explain the reasons for choosing SSP2.6 and SSP8.5 scenarios in this study. Based on CMIP6, there are five shared socioeconomic pathways, and why are the other three scenarios not considered? Although the analysis of all five scenarios is a large workload, as a reader, we need to know what the criteria are for selecting scenarios for this study.\n7. Line 92-93, Line 250-256, Line 333-334: The authors need to clarify the relationship between date of foliar senescence and different droughts, which is very important for selecting the appropriate drought index. Generally, soil drought is not a type of drought, and SSI is generally used to characterize agricultural drought. In meteorological drought and agricultural drought, all kinds of crops on the ground are disaster victims, and hydrological drought is usually expressed as a runoff index. Therefore, is it feasible to select TWSA-DSI when analyzing the correlation between hydrological drought and date of foliar senescence? In my opinion, it would be better to apply SPEI, PDSI, and SRI to represent meteorological droughts, agricultural droughts, and hydrological droughts, respectively."}, "author_response": "Dear Reviewers,\nWe sincerely thank the four Reviewers for your constructive and insightful comments on our\nmanuscript that was originally entitled “Preseason Heatwaves Amplify the Drought Thresholds\nTriggering Earlier Autumn Foliar Senescence over the Northern Hemisphere (ID:\nNCOMMS-24-57320A).” In this revision, we changed the manuscript title as “Drivers of the pre-\nseason drought thresholds triggering earlier autumn foliar senescence in the northern\nhemisphere”. We greatly appreciate the time and effort each of you has dedicated to enhancing the\nquality of our work.\nTo address your valuable feedback, we made comprehensive revisions to the manuscript and\nprovided detailed responses to each of your comments. In our response letter, all the modifications\nand additions are highlighted in blue to clearly indicate how each comment has been addressed.\nSimilarly, corresponding changes in the manuscript are marked in blue, with the referenced line\nnumbers in this letter pertaining to those in the re-uploaded version of the manuscript for your\nconvenience.\nIn our overall response, we summarized the key strategies and major changes implemented to\nimprove the manuscript. This includes expanding the discussions on biophysical mechanisms,\nenhancing the explanations of methodology, and incorporating additional analyses to\ncomprehensively address the concerns of reviewers. In particular, we delved deeper into the\nbiophysical processes that drove our results and clarified our methodological approaches,\nparticularly regarding the use of copula functions with small sample sizes. In addition, we added\nnew analyses to strengthen our findings.\nFor each specific comment, particularly those that required in-depth explanations or\nmethodological adjustments, we provided thorough and detailed responses in the point-by-point\nsection. These detailed explanations aim to clarify the rationale behind our revisions and\ndemonstrate how we have strengthened the robustness and validity of the study. Additionally, we\nincluded new figures and Supplementary Information to support our revised analyses and ensure\nthat our responses are well substantiated.\nWe hope that these revisions effectively address the concerns of the reviewers and enhance\nthe overall quality and clarity of our manuscript. We are committed to further refining our\nmethodology and providing additional validation if required to ensure the reliability and credibility\nof our study."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Thanks a lot for your improvements, however, there are still some minor issues that need to be addressed:\n1. Lines 93-94, there is almost no vegetation in deserts, why are these areas identified PDT-DFS?\n2. In terms of Fig. 1, statistically significant difference analysis should be performed across different biomes and DFS quantiles.\n3. Despite the author incorporate the vegetation hydraulic properties into the analysis, I would also like to know if there is a difference between isohydric and anisohydric biomes, as these two different plant functional types can exhibit different response during drought period. For example, compared to isohydric plants characterized by stomatal closure during drought periods, anisohydric plants can keep their stomata relatively open to maintain photosynthesis, which may exhibit a strong response to drought and heatwave.\n4. Line 104: “The PDT-DFS was not identified for all the regions”, this description seems not correct.\n5. In Fig. S71, almost all trends in the drought thresholds triggering an earlier DFS based on the satellite-derived dataset are not statistically significant, which is quite different from ground observations.\n6. The author highlighted the influence of ecosystem resilience to PDT-DFS. However, the importance of resilience in RF model is very low. In addition, how to calculate the resilience? Most studies used the lag-1 autocorrelation of vegetation optical depth (VOD) to represent the ecosystem resilience, have you used this method to calculate?", "reviewer_2": "I appreciate the efforts authors have made for the extensive revision. Most of y previous suggestions were explained and I have some additional suggestions for authors to consider. After solving these issues, the manuscript could be accepted for publication.\n1. The manuscript needs a very detailed proofreading for clarity. There are many places where grammar mistakes were found and also for descriptions. These should be easy to correct.\n2. The key message should be put forward more clearly. I think the key is the threshold for earlier DFS, not earlier DFS under droughts. Therefore, mechanisms I suggested in the first round is the reason for the thresholds, not for earlier DFS under drought. The latter has already been extensively studied. In this sense, I would like to see some additional analyses, probably because I did not make my points clear. Also, it would be great if the causation analysis could be done to solve these issues that correlation is not always the causation.\n3. The presentation could be improved. Please refer to those figures published for the presentation. For example, the 120 °E could be 90° clock transform.", "reviewer_3": "The revised version addressed my comments well. No further comments.", "reviewer_4": "The application is interesting and the prose in the paper is well-written. All questions in the last version of manuscript have been clarified clearly."}, "author_response": ""}, {"round": 3, "reviewer_comments": {"reviewer_1": "The authors have properly answered my questions. I do not have more comments.", "reviewer_2": "I appreciate the revision authors have done. While although the causality analysis is not involved, the manuscript is promising and could be published after some minor corrections on presentation."}, "author_response": ""}]

Decision letter