Is earlier always better? A comparative assessment of an early-stage response window for rainfall timing in multiyear drought alleviation in China

Multiyear droughts (MYDs) are severe drought events with profound impacts on human activities and ecosystems. However, the timing that maximizes drought attenuation efficiency under additional precipitation (t_optimal) remains insufficiently quantified. With that in mind, we conducted a retrospective grid-based analysis of historical MYDs in China using the Palmer Drought Severity Index (PDSI) over 1961–2020, with a calibration period of 1961–1990, and performed numerical experiments with hypothetical precipitation gradient increases during 1991–2020. We defined k to quantify the sensitivity of drought alleviation to precipitation input. Actual k values derived from the numerical experiments showed strong consistency with theoretical k derived from a drought attenuation quantitative model (DAQM). Theoretical k was further used to quantify MYD attenuation efficiency, with t_optimal identified as the timing of maximum theoretical k (k_max). Overall, drought severity (DS) exhibits a nonlinear response to hypothetical precipitation increases. Both numerical experiment-based and DAQM-derived results identify an early-stage response window within the first three months (t₁∼t₃), and multiple robustness analyses confirm its stability. k_max occurs most frequently at t₁, in 58.79 % of all grids, while t₂ and t₃ were also non-negligible (22.11 % and 11.06 %). This suggests that missing t₁ within this window does not necessarily imply a substantial reduction in k. Compared with humid river basins in southern China, arid and semi-arid northern basins show a higher probability that k at t₂ or t₃ exceed k at t₁. Drought duration (DD) is a key factor, as longer DD is associated with a greater likelihood of t₂ or t₃ being the t_optimal (R²=0.526 and 0.578). We noticed that k is jointly influenced by the cumulative temporal advantage of earlier precipitation input and the higher sensitivity of standardized drought indices under severe drought conditions. This study identifies an early-stage response window during MYD evolution and provides conceptual guidance for understanding how precipitation input timing influences drought alleviation across regions.

1 Introduction

Droughts are recognized as a major natural hazard with profound negative effects on both human activities and the environment (Berdugo et al., 2020; Shi et al., 2021; Palagi et al., 2022). Compared to other natural hazards, droughts often rank higher in key characteristics of losses, such as total loss of life and total economic loss (Mishra and Singh, 2010). For instance, in the United States, drought caused USD 250 billion in damages and nearly 3000 deaths between 1980 and 2020, making it the costliest and second deadliest natural disaster during this period (Ault, 2020; Smith, 2020). In Europe, drought-related financial losses over the three decades prior to 2007 were estimated to exceed EUR 100 billion (Blauhut et al., 2016). In China, prolonged droughts in Southwest China during 2018–2020 resulted in economic losses of USD 240 million, underscoring the intensifying impacts of prolonged droughts (Feng et al., 2025). Such devastating drought events, which typically last more than 12 months (Brunner and Tallaksen, 2019; van Mourik et al., 2025), are referred to as multiyear droughts (MYDs). Under climate change, these persistent, expansive, and difficult-to-recover MYD events have shown an escalating trend in severity, duration, and frequency (Dai, 2013; Zhang et al., 2016; Stevenson et al., 2022; Chen et al., 2025; Wang et al., 2025).

Droughts are primarily related to prolonged reductions in the precipitation for a region (McCabe and Wolock, 2015; Hao et al., 2018; Zhang et al., 2021). For instance, the Millennium Drought in southeast Australia (2001–2009) was the longest drought on record for the region (Van Dijk et al., 2013). Similarly, the California droughts from 2012 to 2015 included the driest three-year period on record in the state (2012–2014) (Mao et al., 2015; Luo et al., 2017). The MYD in Europe from 2014 to 2018 was also record-breaking for the region (Moravec et al., 2021; Büntgen et al., 2021). To better quantify the severity of drought, various drought indices have been developed based on deviations in water availability against long-term normal conditions. Currently, the most widely used drought indices include the Palmer Drought Severity Index (PDSI), Standardized Precipitation Evapotranspiration Index (SPEI), and Standardized Precipitation Index (SPI) (Zhang et al., 2016; Liu et al., 2017; Mukherjee and Mishra, 2021). These indices typically quantify drought severity on a monthly timescale based on precipitation and/or potential evapotranspiration (PET) (Easterling et al., 2007; Vicente-Serrano et al., 2010). They generally standardize monthly effective precipitation time series under the assumption of a normal distribution, enabling the identification of thresholds that delineate different drought intensities (Dai, 2011; Vicente-Serrano et al., 2025). Moreover, by thoughtfully incorporating auto-correlation processes (such as auto-correlation coefficients or the accumulation of antecedent water deficit/surplus), these indices can effectively capture the dynamic evolution of typical drought events, including their onset, development, recovery, and termination (Gupta and Karthikeyan, 2024). To further characterize the dynamic evolution of typical drought events, Yevjevich (1967) proposed the run theory, which extracts drought features based on three key characteristics. These include: (a) Drought duration (DD, the period during which a drought parameter is continuously below the critical threshold), (b) Drought intensity (DI, the maximum deviation of the drought index below the critical threshold during a drought event), (c) Drought severity (DS, the cumulative deficiency of a drought parameter below the critical threshold throughout the event). Such drought event features have been widely applied at station, regional, and global scales to assess drought processes and related impacts (Gu et al., 2020; Wu et al., 2022a; Ullah et al., 2022; Zhou et al., 2023), with DS particularly offering a more integrated measure of drought impact (Cavus and Aksoy, 2020).

Previous studies have examined the dynamic evolution of drought events across the four stages, namely onset, development, recovery, and termination (Ma et al., 2023; Lin et al., 2023). It has been widely observed that the recovery and termination phases are often accompanied by increased precipitation (Haile et al., 2020). However, such hydrometeorological changes typically occur several months after the onset of severe drought, by which time substantial ecological, agricultural, and socioeconomic damages may have already occurred (DeChant and Moradkhani, 2015; Livneh and Hoerling, 2016; Wu et al., 2018). For example, in southeast Australia, the Millennium Drought (2001–2009) was terminated by unusually high rainfall in early 2010, yet this MYD still exhibited prolonged DD and high DS (Yang et al., 2017). Consequently, precipitation occurring earlier in the drought event is more effective in alleviating drought impacts (Fig. 1a and b). From another perspective, drought indices typically represent anomalies from normal conditions in a standardized way, allowing comparisons across space and time. This implies that during extreme drought, when the standardized index deviates far from zero, even small increases in precipitation can induce relatively large numerical shifts in the index, owing to the probabilistic structure of standardized indices. In contrast, when the index is closer to zero, the same precipitation perturbation results in much smaller numerical changes. For instance, using the widely applied PDSI as an example, achieving the same one-unit change in the index requires only about 48 % of the cumulative probability shift under extreme drought conditions (PDSI $=-\mathrm{4}$) compared with that required under milder drought conditions (PDSI $=-\mathrm{1}$) (Fig. 1c and d). This reflects the non-linear probabilistic spacing of the standardized drought indices, whereby values far below the drought threshold are numerically more sensitive to identical perturbations. Consequently, index responses to the same precipitation perturbation can differ substantially depending on the magnitude of PDSI deviations. However, such nonlinear responses are often overlooked. Small increases in precipitation during extreme drought are frequently dismissed as inconsequential, despite their potential to drive meaningful recovery (Pan et al., 2013).

Figure 1 Conceptual and statistical illustrations related to drought alleviation and the normal distribution characteristics of commonly used drought indices (with PDSI as an example). Here, “forcing” refers to a hypothetical incremental water input applied at a given month within a drought event for conceptual illustration. Panels (a) and (b) illustrate the event-level progression of a drought and the associated accumulation of deficit. (a) During the early stage of a drought event, PDSI deviation is relatively small and the accumulated deficit is low. (b) As the drought progresses toward the drought intensity (DI) point, which is defined as the maximum deviation of the index below the critical threshold during the event, the PDSI deviation increases and the accumulated deficit becomes larger. Panels (c) and (d) illustrate the statistical behavior of the standardized index, independent of event-level cumulative deficit. (c) Change in cumulative probability induced by a one-unit shift in PDSI at different levels of PDSI deviation below the drought threshold. (d) PDSI shift required to achieve an equivalent change in cumulative probability at different levels of PDSI deviation.

Currently, the timing that maximizes drought attenuation efficiency under additional precipitation (t_optimal) remains insufficiently quantified. Several existing studies focus on early drought monitoring or alleviation measures during severe drought (Haile et al., 2020). However, we suggest that effective drought alleviation may not necessarily occur at either of these two stages. In particular, cases where early precipitation input and drought conditions are considered together remain under-investigated. Therefore, it is necessary to clarify the relationship between the timing of precipitation input and drought alleviation. Because MYDs involve cumulative water deficits across seasons and years, exhibiting more persistent impacts, the timing of additional precipitation becomes more critical for alleviation, making MYDs particularly suitable for examining precipitation effects on DS. To this end, we propose a new framework involving a series of numerical experiments to retrospectively investigate typical historical MYD events. In this study, MYD events are defined at the grid level as the basic analytical units, with each representing a temporally continuous drought period identified independently within a single grid. Section 2 introduces the study area and datasets. In Sect. 3, representative MYD events are identified based on existing literature, with a particular focus on assessing the extent of drought alleviation during these events. A model based on hypothetical precipitation forcing applied at different months of MYD is developed to parameterize MYD attenuation efficiency. The accuracy of this parameter is evaluated by comparing model-derived estimates with empirical estimates derived directly from the numerical experiments. The results and related discussion are presented in Sects. 4 and 5. We close with a brief conclusion in Sect. 6.

2 Study area and datasets

To assess the role of nonlinear patterns in MYD alleviation, we utilize a daily meteorological dataset suitable for calculating PDSI during the period of 1961 to 2020 over China. This dataset includes precipitation (pr), air temperature (mean, maximum, and minimum), relative humidity, wind speed, and sunshine duration. This dataset consisting of 756 stations is provided by the National Climate Centre of the China Meteorological Administration and has been quality-controlled before being released to the scientific community (https://www.nmic.cn/, last access: 26 August 2025). Given the importance of reliable statistics, we applied additional temporal and spatial consistency control on the data based on the information of the length of available data, the missing data, and the density of stations, etc. First, we conducted temporal consistency control on the dataset. Daily records from 1960–2020 were examined for each station, and stations with more than 5 % missing values were excluded. Short gaps (≤5 d) in the remaining stations were filled using linear interpolation, while longer gaps were replaced with the multi-year mean for the corresponding day. Following this procedure, 620 out of 756 meteorological stations with continuous measurements were selected. These stations are distributed in 10 large river basins covering China, as shown in Fig. 2. Second, we conduct further control on the spatial inhomogeneity of station distribution. We select 199 grid boxes of 2°×2° longitude by latitude (Fig. 2), and each grid box should contain at least one station with continuing meteorological measurements. For the given grid box with multiple stations, the station results were averaged to represent the grid box. In doing that, we reduce the possible influence from temporal and spatial inhomogeneities in the data and ensure that the overall statistics are robust. The resulting grid-level daily meteorological series were then used as the basis for subsequent monthly aggregation and PDSI calculation.

Figure 2 The study area covers 10 large river basins (regions). No. 1, the Songhua River; No. 2, the Liao River; No. 3, Northwestern Rivers; No. 4, the Hai River; No. 5, the Yellow River; No. 6, the Yangtze River; No. 7, the Huai River; No. 8, Southeastern Rivers; No. 9, Southwestern Rivers and No. 10, Pearl River. A total of 199 (2°×2°) grid boxes (color: cyan) were identified as containing at least one station with continuous instrumental meteorological records. In total, 620 meteorological stations (purple dots) were used.

3 Methods

3.1 Overview

To determine t_optimal of drought alleviation, we propose the following methodological framework that evaluates how early-stage hypothetical precipitation forcing and PDSI deviation magnitude jointly modulate the efficiency of alleviating grid-level MYD events. We calculate the PDSI using data from meteorological stations, identify representative grid-level MYD events which capture the key characteristics of persistent and severe droughts, and characterize them in terms of DS, DD, DI, Drought Development Period (DDP), and Drought Intensity Point (DIP, the PDSI value at the drought intensity point) (see Sect. 3.3 for identification criteria). Based on these identified grid-level MYDs, a series of numerical experiments is conducted by applying hypothetical precipitation increases along a prescribed gradient. The numerical experiment consists of the following three steps: (1) Validation of the nonlinear response of DS to incremental precipitation forcing applied along a gradient; (2) construction of a novel drought attenuation quantitative model (DAQM) to describe this response; (3) verification of the applicability of the key parameter derived from the DAQM, namely k, the MYD attenuation efficiency coefficient. k characterizes how rapidly additional precipitation is converted into reductions in DS. For each grid-level MYD, k is evaluated across different forcing timings, and its maximum value (k_max) is identified. The t_optimal is then identified as the timing corresponding to k_max. In addition, we examine the relationship between DD, as a key event characteristic, and MYD attenuation efficiency under hypothetical precipitation forcing. A detailed description of each methodological component is provided in Sect. 3.2–3.4. For clarity, a consolidated table listing principal symbols, acronyms and their description is provided for reference (Table S1 in the Supplement). Figure 3 provides an overview of the workflow and methodological structure of the proposed framework.

Figure 3 Workflow and method overview.

3.2 Estimation of PDSI using the Penman-Monteith method

In this study, PDSI was used to quantify the deviation between observed monthly precipitation and the required precipitation, which was estimated each month under Climatically Appropriate For Existing Conditions (CAFEC). The PDSI was calculated following the CAFEC formulation and calibrated at the grid level over the calibration period 1961–1990 to calculate CAFEC values and establish a fixed climatological reference (Karl, 1986). This ensured that PDSI thresholds remained unchanged throughout the analysis (Fig. S1). The calibration period was intentionally kept independent from the subsequent numerical experiments, which were conducted exclusively for 1991–2020, in order to avoid altering the reference climatology and to allow consistent evaluation of hypothetical precipitation forcing effects (Fig. S1a). During the calibration period, monthly climatological values of the water balance components, including precipitation, PET, recharge, runoff, and losses, were derived to compute CAFEC precipitation, calibration coefficients, and the corresponding Z-index and PDSI time series.

The PDSI formulation was based on a simple two-layer water balance model (Palmer, 1965), which considered both water supply and atmospheric evaporative demand. To ensure the credibility of PDSI results, we used the Penman-FAO56 method with a physically based formulation (Xu et al., 2022) to estimate PET (denoted as PET_pm), as shown in Eq. (1). Specifically, daily PET was first calculated from daily meteorological observations and then aggregated to monthly values for use in the PDSI calculation.

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\mathrm{PET}}_{\mathrm{pm}}=\\ & {\displaystyle \frac{\mathrm{0.408}\cdot \mathrm{\Delta }\cdot \left(R_{\mathrm{n}}-G\right)+\mathit{\gamma }\cdot \frac{\mathrm{900}}{T+\mathrm{273}}\cdot U_{\mathrm{2}}\cdot e_{\mathrm{s}}\cdot \left(\mathrm{1}-\mathrm{Rh}/\mathrm{100}\right)}{\mathrm{\Delta }+\mathit{\gamma }\cdot \left(\mathrm{1}+\mathrm{0.34}\cdot U_{\mathrm{2}}\right)}}\end{array}\end{array}$$

Where Δ is the slope of the vapor pressure curve, R_n is the net radiation (calculated using surface downwelling/upwelling shortwave/longwave radiation), G is the soil heat flux, which was assumed to be negligible (G≈0) at the daily timescale beneath the grass reference surface (Allen et al., 1998), γ is the psychometric constant, T is the monthly averaged air temperature, U₂ is the wind speed at 2 m height, e_s is the saturation vapor pressure at a given air temperature, Rh is the relative humidity.

In this study, long-term PDSI series for the period 1961–2020 were calculated using the tool provided by the Palmer Drought Severity Index (Jacobi et al., 2013).

3.3 Identification and characterization of grid-level MYD events based on run theory

Run theory is a threshold-based method that defines drought events as temporally continuous periods during which PDSI remains below −1 (Yang et al., 2024), such that any month exceeding this threshold marks the termination of an event. According to this method, each identified drought event was further characterized using a set of event-based indicators, including DS, DD, DI, DDP, and DIP.

In this study, the identified grid-level drought events were further selected to capture the key characteristics of MYDs that are both persistent and sufficiently severe (hereafter referred to as typical MYDs), by retaining only those with DD>12 and $\mathrm{DI}<-\mathrm{3}$ (Brunner and Tallaksen, 2019; van Mourik et al., 2025; Chen et al., 2025). MYDs were required to be temporally continuous, with PDSI remaining below −1 throughout the entire event. Since the calibration period was set to 1961–1990 and the CAFEC values needed to remain fixed during numerical experiments, the final selection was restricted to grid-level MYD samples occurring during 1991–2020 (Fig. S1a). A total of 351 grid-level MYD samples meeting these criteria were identified. Each basin contained at least 11 grid-level MYD samples, thereby helping to ensure the statistical robustness of the analysis.

3.4 Numerical experiment design for quantifying drought attenuation efficiency under hypothetical precipitation increase

According to run theory, for a given drought event characterized by its baseline monthly PDSI series (Fig. S2a), a hypothetical precipitation increase applied at a specific month leads to a reduction in DS, assuming all other conditions remain unchanged. This provides a foundation for quantifying drought attenuation efficiency in terms of DS reduction. The numerical experiments therefore apply hypothetical precipitation increases to individual months of a drought event independently, keeping all other months unchanged, and compare the resulting DS responses across months. The timing of hypothetical precipitation increases (t_forcing) is a key experimental parameter, representing the month at which the additional precipitation is applied within a grid-level MYD event.

Since this study focuses on drought alleviation, only hypothetical precipitation increase scenarios with incremental gradients are considered. The procedures of this numerical experiment were given as follows:

Step 1: Validation of nonlinear response

To validate the nonlinear response of MYD alleviation to hypothetical precipitation increases applied at different t_forcing and under varying PDSI values, a single-timing precipitation forcing strategy was adopted. Under this strategy, precipitation is increased at only one t_forcing within an MYD event, while precipitation in all other months remains unchanged. Since the calibration period in this study corresponds to the period 1961–1990, all numerical experiments were conducted for MYDs occurring after 1990, with no temporal overlap between the two periods (Fig. S1a). This design ensures that the reference climate is not influenced by hypothetical precipitation perturbations. Consequently, the baseline remains fixed throughout all gradient-based numerical experiments (Fig. S1a). The forcing strategy for a given MYD event during 1991–2020 is expressed as follows:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\mathrm{pr}}_{\mathrm{forcing}}& ={\mathrm{pr}}_{\mathrm{baseline}}+{\mathrm{pr}}_{\mathrm{rate}}\times \stackrel{\mathrm{‾}}{\mathrm{pr}},{\mathrm{pr}}_{\mathrm{rate}}\\ & =\mathrm{0}\mathit{\%},\mathrm{5}\mathit{\%},\mathrm{10}\mathit{\%},\mathrm{\dots },\mathrm{100}\mathit{\%}\end{array}\end{array}$$

Where pr_forcing is the precipitation in the targeted month after applying the hypothetical forcing, pr_baseline is the observed precipitation in the same month under the baseline (unforced) scenario,$\stackrel{\mathrm{‾}}{\mathrm{pr}}$ is the long-term mean monthly precipitation, ${\mathrm{pr}}_{\mathrm{rate}}\left(=\frac{{\mathrm{pr}}_{\mathrm{forcing}}-{\mathrm{pr}}_{\mathrm{baseline}}}{\stackrel{\mathrm{‾}}{\mathrm{pr}}}\right)$ is the relative change in precipitation, and ${\mathrm{pr}}_{\mathrm{rate}}\times \stackrel{\mathrm{‾}}{\mathrm{pr}}$ is the imposed hypothetical precipitation increment. The long-term mean monthly precipitation was used as the reference quantity for defining the perturbation magnitude because it provides a physically interpretable measure of local climatological conditions and is highly consistent with the standard deviation of monthly precipitation across the study area (Fig. S3). Considering that MYD responses may vary under different precipitation perturbation magnitudes, pr_rate ranging from 0 % to 100 % at 5 % increments were used to systematically evaluate MYD responses across different perturbation conditions, thereby integrating both small and relatively large hypothetical precipitation perturbations. The upper bound of 100 % was selected to ensure that, even under severe antecedent drought conditions, the additional precipitation could bring monthly precipitation to or above its long-term mean and potentially shift conditions toward a relatively wet state.

As described above, applying precipitation forcing alters the DS of the corresponding MYD event. Therefore, the relationship between pr_forcing and the resulting change in DS provides a basis for examining the nonlinear response of drought alleviation to hypothetical precipitation increases.

Step 2: Construction of DAQM

To further parameterize drought attenuation efficiency in a consistent and comparable way across grids, forcing timings and MYDs, we developed DAQM. DAQM describes the relationship between hypothetical precipitation perturbations and the resulting reductions in event-level DS. It provides a unified framework for estimating the attenuation efficiency coefficient (k) while accounting for the nonlinear responses implied by Conceptual Models 1 and 2.

DS is defined as the cumulative sum of PDSI values over DD of a grid-level MYD period. In the model formulation, attenuation of DS under hypothetical precipitation increases is represented as a nonlinear diminishing-return process. This functional form is motivated by numerical experiments, which consistently show that incremental precipitation inputs produce progressively smaller marginal reductions in DS across grids and forcing timings, indicating a bounded drought alleviation response.

To represent this behavior, we adopt a first-order kinetic-like formulation, in which the instantaneous rate of DS reduction is proportional to the remaining DS. This structure captures rapid initial alleviation when residual DS is high and progressively slower alleviation as DS approaches zero.

Accordingly, Eq. (3) expresses the drought attenuation process in differential form, linking the instantaneous rate of DS reduction to the residual DS:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial {\mathrm{DS}}_{\mathrm{rate}}}{\partial {\mathrm{pr}}_{\mathrm{rate}}}}=k\cdot \left(\mathrm{1}-{\mathrm{DS}}_{\mathrm{rate}}\right)\end{array}$$

Where ${\mathrm{DS}}_{\mathrm{rate}}(=\frac{{\mathrm{DS}}_{\mathrm{baseline}}-{\mathrm{DS}}_{\mathrm{forcing}}}{{\mathrm{DS}}_{\mathrm{baseline}}}$, DS_baseline is the observed DS under the baseline scenario, and DS_forcing is the DS after the hypothetical precipitation is applied) is the relative change in DS.

k is the attenuation efficiency coefficient, representing the proportionality with which additional precipitation is converted into reductions in residual DS. Because k is dimensionless, it allows direct comparison of attenuation efficiency across different t_forcing and regions.

It is clear that when no hypothetical precipitation forcing is applied (pr_rate=0), there is no change in DS (DS_rate=0). Likewise, when sufficient forcing is applied, the MYD is fully alleviated (DS_rate=1). Thus, the model is subject to the following boundary and limiting conditions:

$$\begin{array}{}\text{(4)}& {\displaystyle }\left\{\begin{array}{r}{\displaystyle }{\mathrm{DS}}_{\mathrm{rate}}\left(\mathrm{0}\right)=\mathrm{0}\\ {\displaystyle }\underset{{\mathrm{pr}}_{\mathrm{rate}}\to \mathrm{\infty }}{lim}{\mathrm{DS}}_{\mathrm{rate}}\left({\mathrm{pr}}_{\mathrm{rate}}\right)=\mathrm{1}\end{array}\right.\end{array}$$

Thus, based on the differential Eq. (3) above and the boundary and limiting conditions in Eq. (4), the analytical solution describing the relationship between DS_rate and pr_rate can be obtained as follows:

$$\begin{array}{}\text{(5)}& {\mathrm{DS}}_{\mathrm{rate}}=\mathrm{1}-\exp \left(-k\cdot {\mathrm{pr}}_{\mathrm{rate}}\right)\end{array}$$

To ensure that this exponential structure reflects observed system behavior rather than an assumed mathematical form, we compared formulations inspired by zero-order, first-order, and second-order kinetic processes using historical MYDs during 1991–2020. Across most grids and events, the first-order structure more accurately reproduced the curvature observed in the numerical experiments (Table S2 in the Supplement). Extension to an n-order formulation showed that the estimated n converged toward unity (n≈0.98), further supporting the first-order representation. We also compared the exponential form against alternative empirical relationships. Logistic curves did not match the observed pattern. Power-law functions were unable to satisfy the required boundary conditions.

Step 3: Verification of parameter robustness

To evaluate the robustness of the drought attenuation efficiency parameter k, we compared model-derived estimates obtained from DAQM (k_theoretical) with empirical estimates calculated directly from the numerical experiments (k_actual). Specifically, k_theoretical is obtained by fitting the DAQM analytical solution to the full set of experimental outcomes, thereby capturing the nonlinear and bounded response of DS to hypothetical precipitation increases. In contrast, k_actual represents the empirical sensitivity of DS to small precipitation perturbations and is approximated from the initial slope of the relationship between DS_rate and pr_rate under the smallest forcing step ($k_{\mathrm{actual}}=\frac{{\mathrm{DS}}_{\mathrm{rate}}\left({\mathrm{pr}}_{\mathrm{rate}}=\mathrm{0.05}\right)}{\mathrm{0.05}}$, where 0.05 is the smallest perturbation step used in the experiments).

The relative error between k_theoretical and k_actual was used as a diagnostic measure to assess the numerical reliability of the model-based parameterization:

$$\begin{array}{}\text{(6)}& \mathit{ϵ}=\left({\displaystyle \frac{k_{\mathrm{theoretical}}-k_{\mathrm{actual}}}{k_{\mathrm{actual}}}}\right)\times \mathrm{100}\mathit{\%}\end{array}$$

Where ϵ is the relative error, k_theoretical is the model-derived attenuation efficiency coefficient, and k_actual is the empirical estimate obtained from numerical experiments.

In subsequent analyses, k_theoretical is used as the primary indicator because it captures the overall nonlinear response structure of DS to hypothetical precipitation increases and provides a more comprehensive and comparable measure of drought attenuation efficiency across grids and forcing timings. k_actual is employed solely to verify internal consistency and is not used in other analyses.

4 Results

4.1 Response relationship between DS and hypothetical precipitation increase in China

Focusing on drought alleviation, the numerical experiments were constrained to the t_forcing range of DDP, with the minimum DDP among the selected MYD events being 8 months. Accordingly, t₁, t₂, …, and t₈ represent the first to the eighth month of MYD. For each t_forcing from t₁ to t₈, DS_rate values were calculated across hypothetical precipitation perturbations with pr_rate ranging from 0 to 1 at 0.05 increments. Figure 4a presents the empirical DS_rate responses derived directly from the numerical experiments, summarized as the mean ± standard error (SE) across MYD events under different precipitation perturbation scenarios, mainly to illustrate the overall attenuation behavior across the study area. The resulting responses were then fitted for each t_forcing within each MYD event using the DAQM analytical form described in Sect. 3.4 (Eq. 5). The fitted relationships closely reproduced the numerical experiment outputs at the event scale with R² values generally exceeding 0.90, confirming the suitability of the DAQM analytical form. This fitting yielded one theoretical attenuation efficiency parameter (k) for each t_forcing, representing the drought attenuation efficiency of a given event at a specific timing. The corresponding k values are denoted as k₁, k₂, …, and k₈, respectively.

Figure 4 DS responses to hypothetical precipitation perturbations and cumulative absolute PDSI from event onset to a given month within a single event relative to the DS of that event during DDP in China. (a) Mean responses under different hypothetical precipitation perturbations for t₁∼t₈ across all MYD events, with error bars representing standard error. (b) Temporal evolution of $\frac{\mathrm{cumulative}\left|\mathrm{PDSI}\right|}{\mathrm{DS}}$ during DDP across all MYD events.

To preserve regional differences, MYD events were identified and analyzed separately at each grid across the study area, rather than being calculated from regionally averaged conditions. By retaining individual grid-level responses, the analysis preserves spatial heterogeneity while comprehensively reflecting the overall MYD response characteristics across the study area. Based on these grid-level events, the changes in DS were not proportional along the increasing precipitation gradient. Within the same pr_rate range from 0 to 1, all timings from t₁ to t₈ exhibited a similar nonlinear diminishing-return response pattern (Fig. 4a). The overall responses for t₁, t₂, and t₃ remained relatively close, whereas a clearer reduction emerged from t₄ onward. Effect sizes between two groups were evaluated using Cohen's d, and the magnitude of differences was quantified using its absolute value (hereafter referred to as Cohen's d). Pairwise statistical comparisons further showed that although the differences among t₁∼t₃ were statistically significant (p<0.001), the corresponding effect sizes were extremely small (Cohen's $d=\mathrm{0.002}\sim \mathrm{0.003}$) (Table S3). In contrast, comparisons involving t₄ exhibited substantially larger effect sizes (Cohen's $d=\mathrm{0.231}\sim \mathrm{0.342}$) (Table S3). These results indicate that, when aggregated across all MYD events over the study area, the overall drought attenuation efficiencies associated with t₁∼t₃ remained relatively close before a clearer decline emerges from t₄ onward.

Since hypothetical precipitation inputs applied at a given t_forcing can only influence the subsequent stage of MYD evolution, cumulative deficits were further calculated to quantify the drought impacts accumulated before each timing. Because PDSI values during MYD events remain below the drought threshold ($\mathrm{PDSI}<-\mathrm{1}$), cumulative absolute PDSI was used in the analysis. Based on this, the ratio of cumulative absolute PDSI from event onset to a given month within a single event relative to the DS of that event $\left(\frac{\mathrm{cumulative}\left|\mathrm{PDSI}\right|}{\mathrm{DS}}\right)$ was defined to quantify the proportion of drought impacts accumulated before each t_forcing. The ratio increased progressively with MYD development across all MYD events (Fig. 4b). The $\frac{\mathrm{cumulative}\left|\mathrm{PDSI}\right|}{\mathrm{DS}}$ values corresponding to t₁, t₂, t₃, and t₄ were approximately 0.04, 0.10, 0.17, 0.25, respectively. These results indicate that a progressively larger fraction of DS had already accumulated by later t_forcing during MYD development, thereby reducing the remaining alleviation potential.

Within the relatively early MYD stages (t₁∼t₃), further basin-scale analyses of pairwise comparisons among t₁∼t₃ showed that the distinctions among early t_forcing became substantially more pronounced (Fig. S4). The median Cohen's d values for pairwise comparisons among t₁∼t₃ generally exceeded 0.80 in all ten river basins, corresponding to large effect sizes, while the minimum values were mostly greater than 0.50, indicating at least medium effect sizes. These results indicate that the differences among t₁∼t₃ remain distinguishable at the river basin scale, supporting further investigation of attenuation efficiency differences among early t_forcing.

4.2 Validation of the MYD attenuation efficiency coefficient in China

For each grid, the t_optimal was first identified based on the DAQM-derived k. Figure 5a shows the spatial distribution of the corresponding t_optimal across the study area. Overall, t₁ dominates, accounting for 58.79 % of the grids, whereas t₂ and t₃ account for 22.11 % and 11.06 %, respectively. To further evaluate the robustness of the DAQM-derived t_optimal, an independent estimate was obtained from DS_rate responses (Fig. S5). During the numerical experiments, it was found that the identified t_optimal for a given MYD event was not always identical across different hypothetical precipitation perturbation scenarios, indicating that precipitation perturbation magnitude itself can influence the identified t_optimal. To reduce the influence of perturbation magnitude on t_optimal identification, 21 hypothetical precipitation perturbation scenarios (pr_rate=0, 0.05, 0.10, …, 1.00) were separately evaluated for each MYD event, thereby integrating both small and relatively large perturbations into a comprehensive evaluation method. For each scenario, the corresponding optimal timing was first determined according to the maximum DS_rate among different t_forcing. The mode of the resulting 20 optimal timings was then used to represent t_optimal across different precipitation perturbation magnitudes (within the mode-based identification process, scenario without additional precipitation input, i.e., pr_rate=0, was not included because identical DS_rate across all t_forcing prevented a unique optimal timing from being identified), hereafter referred to as the mode-based approach. Because t_optimal is expressed as a discrete positive integer, the mode was adopted instead of the mean or median to ensure that the identified t_optimal remained physically interpretable as an integer month. Notably, the DAQM-based estimates are also derived from all hypothetical precipitation perturbation scenarios through least-squares fitting, making the mode-based and DAQM-derived results directly comparable. The mode-based approach yields proportions of t₁ (67.84 %), t₂ (12.06 %), and t₃ (14.57 %) (Fig. S5). Comparison between the two approaches shows that 75.88 % of the grids exhibit identical t_optimal (Fig. 6b). Further analyses focusing only on grids where the two approaches yield identical t_optimal show proportions of t₁ (72.84 %), t₂ (11.92 %), and t₃ (10.60 %) (Fig. S6).

Figure 5 Validation and robustness assessment of the MYD attenuation efficiency coefficient (k). (a) Spatial distribution of the DAQM-derived t_optimal for each grid. The inset shows the frequency histogram of the DAQM-derived t_optimal across all grids. (b) Spatial distribution of t_optimal differences (mode-based t_optimal minus DAQM-derived t_optimal). The inset shows the frequency histogram of timing differences between mode-based and DAQM-derived t_optimal across all grids. (c) Spatial distribution of ϵ for each grid. (d) Spatial distribution of t_optimal differences (small perturbation-based t_optimal minus DAQM-derived t_optimal). The inset shows the frequency histogram of timing differences between small perturbation-based and DAQM-derived t_optimal across all grids. (e–g) Median actual versus median theoretical k values for the ten major river basins in (e) t₁, (f) t₂, and (g) t₃. The red dashed line is the 1:1 line. The black solid line is the fitted line based on values from 10 river basins. Orange crosses are the values for individual river basins.

Figure 6 Spatial distribution of the changes in k (Δk) derived from the DAQM model. (a) Change in k between t₂ and t₁ (Δk₂₁). (b) Change in k between t₃ and t₁ (Δk₃₁).

According to the DAQM constructed in Sect. 3.4, as the increment in pr_rate approaches zero, the corresponding DS_rate also approaches zero. Under this condition, the estimated value of k approaches the actual k. Therefore, the actual values of k can be approximated using the numerical experiment controlling a small perturbation of a individual input variable (pr_rate=0.05) and its corresponding DS_rate. Based on this principle, actual values of k were calculated from the numerical experiments and compared with the theoretical k values derived from the DAQM developed in Sect. 3.4, at the grid scale. Figure 5c shows the spatial distribution of ϵ across all given grids, with 72.86 % of the grids exhibiting ϵ values below 5 %. Slightly larger deviations mainly occur in parts of northern and southern China, whereas semi-arid regions generally exhibit relatively stronger agreement between theoretical and actual estimates. Despite these regional differences, the overall ϵ remains modest, with most grids exhibiting values around 6 %, confirming the robustness of the DAQM across the study area (Figs. 5c and S4). To further evaluate whether ϵ influences t_optimal, the small perturbation experiments (pr_rate=0.05) were also used to independently identify t_optimal (Fig. S7). Comparison with the DAQM-derived t_optimal shows agreement in 89.95 % of the grids (Fig. 5d), indicating that the ϵ does not substantially alter the spatial pattern of t_optimal. To further evaluate the robustness of the DAQM framework at the basin scale, comparisons between theoretical (DAQM-derived) and actual k values (small perturbation-based) were conducted across the ten major river basins in China. A strong positive correlation was observed between the theoretical and actual values, with R² values exceeding 0.90 for t₁, t₂, and t₃ across the ten river basins in China (Fig. 5e–g). The comparison results for each month revealed no significant deviation from the 1:1 reference line, with slopes of 1.063, 1.060, and 1.037 for t₁, t₂, and t₃, respectively (Fig. 5e–g).

Overall, these results demonstrate the robustness and applicability of the theoretical values at both the grid and river basin scales. For this reason, the DAQM was further employed to investigate drought attenuation efficiency at different t_forcing. This quantitative model was also used to identify t_optimal for precipitation input across t₁, t₂, and t₃.

4.3 Quantification of the relationship between MYD attenuation efficiency and precipitation input timing in China

To identify t_optimal in various river basins during t₁, t₂, and t₃, a series of numerical experiments were conducted using the DAQM. Figure 6 presents the spatial distributions of the changes in k between t₂ and t₁ ($\mathrm{\Delta }{k}_{\mathrm{21}}=k_{\mathrm{2}}-k_{\mathrm{1}}$), and between t₃ and t₁ ($\mathrm{\Delta }{k}_{\mathrm{31}}=k_{\mathrm{3}}-k_{\mathrm{1}}$). As demonstrated by the validation results shown in Sect. 4.2, the DAQM framework performs robustly across most regions, supporting the interpretation of the following spatial patterns. Additional analyses based on small perturbation experiments (pr_rate=0.05) further confirm that the DAQM-derived results generally capture the spatial patterns (Fig. S8). Strong consistency was observed between the responses derived from the small perturbation experiments and those obtained from the full hypothetical precipitation-gradient framework, with correlation coefficients of R=0.957 for Δk₂₁ (Figs. 6a and S8a) and R=0.955 for Δk₃₁ (Figs. 6b and S8b) (p<0.001 in both cases). Therefore, the spatial patterns shown in Fig. 6 can be considered reliable for interpreting drought attenuation efficiency.

To enable a clearer comparison of magnitude among cases where Δk<0, absolute values were used. Notably, the larger absolute values of negative Δk (Δk<0) were mainly observed in the humid river basins of southern China, particularly in the Pearl River and Yangtze River basins (Fig. 6 and Table 1). By comparison, smaller absolute values of negative Δk were detected in the arid and semi-arid northern regions of China, particularly in the Liao River basin (Fig. 6 and Table 1). Moreover, in 118 out of 199 grids (59.30 %), the absolute values of negative Δk₃₁ were larger than those of Δk₂₁ (Fig. 6 and Table 1). This disparity indicates that, in more than half of MYD occurrences, the greatest attenuation efficiency tends to occur at earliest t_forcing (t₁). Interestingly, we found that 40 out of 199 grids (20.10 %) exhibited positive values in both Δk₂₁ and Δk₃₁ (Fig. 6). These results indicate that, in a non-negligible portion of grids, the greatest attenuation efficiency does not occur at t_1.

Table 1 Summary of the number of grids, median values of DD, and changes in theoretical k (derived from the DAQM model; $\mathrm{\Delta }{k}_{\mathrm{21}}=k_{\mathrm{2}}-k_{\mathrm{1}}$, $\mathrm{\Delta }{k}_{\mathrm{31}}=k_{\mathrm{3}}-k_{\mathrm{1}}$, and $\mathrm{\Delta }{k}_{\mathrm{32}}=k_{\mathrm{3}}-k_{\mathrm{2}}$) and actual k (derived from the small perturbation experiments with pr_rate=0.05), across the ten river basins and all of China.

Note: The bold value marks the maximum P₅₀ among the ten river basins in China, and the bold-italic value marks the minimum P₅₀. P₅₀ values are derived from both the DAQM-derived and small perturbation-based (pr_rate=0.05) results.

These experiments also aimed to identify the key drought event characteristics responsible for the variation in t_forcing. Since DS represents the cumulative water deficit over the entire DD, and k is derived from DS, DD constitutes an important factor influencing k in addition to the standardized characteristics of PDSI itself. Therefore, the spatial distribution and statistical characteristics of DD are examined to help interpret basin-scale differences in k (Fig. 7 and Table 1). We summarize the probability density function (PDF) of Δk₂₁, and Δk₃₁ for all of grid-level MYD events across 10 river basins in China (Fig. 7 and Table 1). The Δk values exhibit a roughly normal distribution across river basins (Fig. 7). The median values of both Δk₂₁ and Δk₃₁ are negative in 9 out of 10 basins, with the median of Δk₃₁ being lower than that of Δk₂₁ (Table 1). In most basins, the standard deviation of Δk₃₁ is also larger than that of Δk₂₁, indicating greater variability in the efficiency difference between t₃ and t₁ (Fig. 7). These results suggest an efficiency advantage of early precipitation input (t₁>t₂ and t₁>t₃). This advantage becomes more pronounced when comparing t₁ and t₃. However, positive Δk₂₁ and Δk₃₁ values were found in 33.67 % and 23.12 % of all grids, respectively (Fig. 6a and b). Similar results were also evident in the PDFs for all of China and each individual basin (Fig. 7). Consistent patterns are observed in both the DAQM-derived results (Fig. 7) and the small perturbation-based results (Fig. S9). These results suggest a non-negligible role of the higher sensitivity of normally distributed standardized drought indices to identical precipitation perturbations under more severe drought conditions. Longer DD tends to be associated with Δk values closer to or even above zero (Fig. 7 and Table 1). In particular, the Liao River, with a relatively long DD of 28 months (DD_{all China}=21 months), shows median Δk values of 0.022 for Δk₂₁ and 0.003 for Δk₃₁ (Fig. 7 and Table 1). In contrast, the Pearl River, characterized by a shorter DD of 19 months, exhibits more negative median Δk values of −0.162 and −0.200 for Δk₂₁ and Δk₃₁ (Fig. 7 and Table 1), respectively.

Figure 7 Spatial distribution of the median drought duration from 1991 to 2020 across the 10 river basins in China. Surrounding panels show the probability distribution of Δk derived from the DAQM model for all of China and each basin, based on grid-level MYD samples. Specifically, the orange lines represent Δk₂₁, and the green lines represent Δk_31. The x-axis in each panel denotes Δk, while the y-axis denotes probability. The number in the upper-right corner of each panel indicates the basin code.

To further investigate the influence of DD as an important event characteristic, an equal-width binning method was applied to examine the variation of Δk under different DD values. Specifically, the original DD values of the selected MYD events ranged from 13 to 51 months. These events were grouped into bins according to their DD values. To ensure the robustness of the results (each bin contained no fewer than 10 grid-level events), the DD range was limited to 13 to 40 months, and the bin width was set to 4 months. This setting yielded a total of seven bins (Fig. 8). We calculated the median values of both Δk and DD within each bin. As shown in Fig. 8, we found a positive correlation between DD and Δk, with R² values of 0.526 and 0.578 for Δk₂₁ and Δk₃₁, respectively. Notably, when DD exceeds 29 months (as shown in Fig. 8), the Δk values approaching zero, indicating that the drought attenuation efficiencies of t₂ and t₃ become comparable to that of t₁. In contrast, when DD is shorter (around 17 months, as shown in Fig. 8), t₁ demonstrates notably higher drought attenuation efficiency than t₂ and t₃. These results suggest that DD is a key factor influencing t_optimal.

Figure 8 Median Δk of MYD under different DD. The orange line represents Δk₂₁, and the green line represents Δk_31. Shaded areas indicate half the standard error as the confidence interval.

To clarify the alleviation advantages of the t_optimal, we consider both the cumulative temporal advantage of earlier precipitation input and the statistical properties of standardized drought indices. Figure 9 and Table 2 present a comparison between the drought attenuation efficiency of t_optimal and those of t₁ and t_DIP (the month corresponding to DIP within a MYD event). The corresponding k at t₁ (k₁), t_DIP (k_DIP), and t_optimal (k_max) exhibit differences. For simplicity, two ratios, $k_{max}/k_{\mathrm{1}}$ and $k_{max}/k_{\mathrm{DIP}}$, are used to quantify the differences in k between t_optimal and t₁, and between t_optimal and t_DIP, respectively. In terms of magnitude, $k_{max}/k_{\mathrm{1}}$ is typically close to 1, while $k_{max}/k_{\mathrm{DIP}}$ consistently exceeds 2 (Table 2). From the perspective of the cumulative temporal advantage of earlier precipitation input, t₁ tends to dominate (Fig. 9a and Table 2). Nevertheless, the presence of grids with $k_{max}/k_{\mathrm{1}}>\mathrm{1}$ suggests t_optimal>1 should not be overlooked (Fig. 9a). At the basin scale, t₁dominates in most regions, except in the Liao River basin, where the median $k_{max}/k_{\mathrm{1}}$ reaches 1.135 (Table 2). In contrast, t_optimal offers a more substantial advantage in drought attenuation efficiency compared to t_DIP (Fig. 9b and Table 2). This advantage is particularly pronounced in the Liao River and Songhua River basins, where the median $k_{max}/k_{\mathrm{DIP}}$ values reach 15.430 and 12.555, respectively. By contrast, the advantage of t_optimal is less evident in the Southeastern River and Hai River basins, with median $k_{max}/k_{\mathrm{DIP}}$ values of 3.428 and 2.413, respectively.

Figure 9 Spatial distribution of the ratios of k derived from the Optimal Solution Model relative to those derived separately from Conceptual Model 1 and Conceptual Model 2. (a) The ratio of k between the t_optimal and t₁ ($k_{max}/k_{\mathrm{1}}$), (b) the ratio of k between the t_optimal and t_DIP ($k_{max}/k_{\mathrm{DIP}}$).

Table 2 Summary of the median values of $k_{max}/k_{\mathrm{1}}$ and $k_{max}/k_{\mathrm{DIP}}$ across the ten river basins and all of China.

Note: The bold value marks the maximum P₅₀ among the ten river basins in China, and the bold-italic value marks the minimum P₅₀.

5 Discussion

5.1 MYD alleviation response to hypothetical increased precipitation

Current studies on precipitation-based drought alleviation, referring to the dominant role of precipitation in driving drought alleviation, have primarily focused on the early stage of drought development (DeChant and Moradkhani, 2015; Livneh and Hoerling, 2016; Wu et al., 2018). Yang et al. (2017) integrated early precipitation input with drought-induced losses and concluded that earlier precipitation is more effective in alleviating drought impacts. This interpretation is consistent with the cumulative temporal advantage associated with earlier precipitation input considered in this study. However, due to the statistical properties of standardized drought indices, more severe drought conditions are often associated with higher attenuation efficiency under identical precipitation perturbations (Fig. 1). To further interpret the mechanisms controlling drought attenuation efficiency, two conceptual models were considered (Fig. S2). Conceptual Model 1 considers the same hypothetical precipitation increase applied at different t_forcing within an MYD event (Fig. S2b). Because the precipitation perturbation can influence DS only from t_forcing onward, earlier t_forcing allows a longer period during which the drought process can be affected. Consequently, t_optimal under Conceptual Model 1 tends to occur at the beginning of the event (t₁), i.e., t_optimal→t₁. Conceptual Model 2 features the same precipitation forcing applied at the same t_forcing, but under different PDSI values (Fig. S2c). As a standardized index, PDSI (similar to widely used normally distributed standardized indices such as SPEI and SPI) exhibits tail sensitivity, whereby values further below the drought threshold tend to experience larger numerical responses under identical precipitation perturbations. Because DS is calculated as the cumulative sum of monthly PDSI values over DD, this property of PDSI also affects event-level DS. This mechanism forms the basis of Conceptual Model 2 and favors t_optimal toward the month of DIP occurrence (t_DIP), i.e., t_optimal→t_DIP. The interaction between these two mechanisms suggests that t_optimal should theoretically fall within the interval from t₁ to t_DIP, because t_DIP of MYD events typically occurs much later than the onset (t_DIP≫t₁) (Fig. S2d). The nonlinear pattern in DS response to increased precipitation (Figs. 1 and 4) can be more comprehensively understood by jointly considering these two conceptual models. To further clarify this relationship, attenuation efficiencies at different t_forcing were quantified using the DAQM developed in Sect. 3.4. Under the same hypothetical precipitation forcing, a higher value of k at a given t_forcing indicates greater drought attenuation efficiency. The results further show that attenuation efficiencies among t₁∼t₃ remain relatively close, whereas clearer separation emerges after t₄ (Fig. 4a and Table. S3). This pattern suggests that t₁∼t₃ constitute an important early-stage response window during which the effects described in Conceptual Models 1 and 2 partially offset each other. However, progressively increasing $\frac{\mathrm{cumulative}\left|\mathrm{PDSI}\right|}{\mathrm{DS}}$ during drought development (e.g., approximately 0.04, 0.10, 0.17, and 0.25 at t₁, t₂, t₃, and t₄, respectively) indicate the continuous accumulation of drought impacts during MYD evolution. Once $\frac{\mathrm{cumulative}\left|\mathrm{PDSI}\right|}{\mathrm{DS}}$ reaches a sufficiently high level (approximately 0.25 at t₄), the cumulative temporal advantage associated with earlier precipitation input becomes dominant. These findings further reflect the complexity of drought alleviation processes and can be interpreted through the interaction between the two conceptual models.

Consistent with several existing studies (Pan et al., 2013; Haile et al., 2020), t₁ remains the dominant t_optimal across China. However, we noted that t₂ and t₃, while occurring later than t₁, are still non-negligible in the early stage of MYD. Spatially, the grids where t₂ and t₃ are dominant account for 22.11 % and 11.06 % of all grids, respectively (Fig. 5a). According to the classification of climate zones based on average annual precipitation (Yang et al., 2022), compared with the humid river basins in southern China, the arid and semi-arid northern regions exhibited a higher probability that drought attenuation efficiency at t₂ and t₃ exceeds that at t₁ (Figs. 6 and 9a). This interpretation is supported by the strong agreement between theoretical and actual k values shown in Fig. 5. We further noticed that the Liao River basin (semi-arid basin) and Pearl River basin (humid basin) exhibited relatively high and low MYD attenuation efficiency, respectively. This difference may be associated with the generally more severe drought conditions in the Liao River basin compared to those in the Pearl River basin (Huang et al., 2019; Lai et al., 2024; Zhang et al., 2025). Due to the statistical properties of standardized drought indices, increased precipitation tends to induce stronger attenuation responses under more severe drought conditions. Consequently, the influence of the second mechanism becomes relatively more evident in arid and semi-arid regions, increasing the probability that t₂ or t₃ may exceed t₁ in attenuation efficiency within the early-stage response window. By contrast, in humid basins characterized by relatively milder drought conditions, the cumulative temporal advantage associated with earlier precipitation input remains dominant, resulting in a stronger preference for t₁. This conclusion was consistent with previous findings (Wang et al., 2025), showing the importance of earlier precipitation input for alleviating MYD in the humid basins of China.

Overall, our results do not suggest that “later is better” for MYD alleviation. Instead, they reveal the existence of an early-stage response window during which t₂ and t₃ may still exhibit attenuation efficiencies comparable to or even exceeding that of t₁. This finding is particularly important for MYDs because these events involve cumulative water deficits across seasons and years, making the timing of precipitation input especially critical for drought alleviation. Missing t₁ does not necessarily imply that the opportunity for maximum attenuation efficiency has already been lost, provided that precipitation input still occurs within the early-stage window.

5.2 Robustness of optimal timing identification within the early-stage response window

The basin-scale Cohen's d results indicate that the attenuation efficiency differences among t₁∼t₃ are distinguishable and cannot be attributed to random variability or noise (Fig. S4). Accordingly, the robustness of the DAQM-derived t_optimal was evaluated through comparisons among the DAQM-derived results, the mode-based results, and the small perturbation-based results. Despite methodological differences among these three approaches, all consistently indicate that t₁ dominates overall, while the possibility that t₂ and t₃ may exceed t₁ within the early-stage response window remains non-negligible (Figs. 5a, S5, S6, and S7).

Comparison between the mode-based and DAQM-derived approaches shows that 75.88 % of grids exhibit identical t_optimal (Fig. 5b). The distribution of timing differences is approximately symmetric and centered around zero, suggesting that the discrepancies are spatially random rather than systematically biased. In addition, Fig. S6 shows that, even when only grids with identical results are considered, t₂ and t₃ still account for a substantial proportion of t_optimal (11.92 % and 10.60 %, respectively), confirming their non-negligible role within the early-stage response period. Generally, small perturbation experiments (pr_rate=0.05) represent more frequently occurring weak perturbation conditions in real-world situations. The comparison between small perturbation-based and DAQM-derived estimates further supports the robustness of the identified t_optimal. Although relative errors in k exist across grids, the overall relative error magnitude remains small across China, with 72.86 % of grids showing relative errors in k below 5 % and most of the remaining grids around 6 %, indicating relatively limited deviations overall. Even though 27.14 % of grids exhibit relative errors exceeding 5 %, this has limited influence on the identification of t_optimal, as 89.95 % of grids show identical results between the two approaches (Fig. 5d). The distribution of t_optimal differences in Fig. 5d further indicates that the remaining discrepancies are not systematically biased. In addition, compared with the small perturbation-based results, the DAQM-derived probability distributions of Δk appear smoother (Figs. 7 and S9), reflecting that the DAQM integrates attenuation responses across the full precipitation-gradient experiment.

Overall, the DAQM-derived results demonstrate strong robustness across multiple validation approaches. The DAQM provides a theoretical basis that helps interpret and justify the patterns observed in the numerical experiments, thereby offering a more physically interpretable basis for attenuation efficiency and t_optimal identification.

5.3 Key drought event characteristics influencing MYD attenuation efficiency

Recent studies have shown that droughts are increasing in both frequency and severity (Dai et al., 2020; Fischer et al., 2025), although substantial regional differences exist. Similarly, Chen et al. (2025) highlighted that MYDs are also becoming more frequent and severe at global scales. To further characterize the dynamic evolution of MYD, three commonly used features include DD, DS, and DI. Among these, DD is a key factor for identifying MYDs, as MYDs are typically defined by a duration exceeding 12 months (Massari et al., 2022). DS represents the overall impact of MYD, while DI is related to the onset timing of the recovery stage. In current studies, standardized drought indices are commonly used to quantify drought severity (Mishra and Singh, 2010). According to the statistical characteristics of the normal distribution (Dai, 2011), drought indices responses to identical precipitation inputs can vary considerably depending on the magnitude of index values below the drought threshold. This aligns with Conceptual Model 2 presented in Sect. 5.1, where small increases in precipitation have the potential to induce substantial shifts toward alleviation when the index deviates far below the drought threshold. Our findings further indicate that t₁ and t_DIP of MYD do not coincide, especially in longer-duration MYD, where t_DIP tends to occur much later than t₁. With that in mind, we examined the relationship between DD and k. Our results demonstrate a positive correlation between DD and k in MYD (Figs. 7 and 8). Specifically, basins in northern China, characterized by longer DD, tend to have a higher probability of t₂ or t₃ being the t_optimal within the early-stage response window, in contrast to southern basins with shorter DD (Fig. 7 and Table 1).

5.4 Limitations and uncertainties

Droughts are not only influenced by pr, but also by other non-precipitation factors such as wind speed, air temperature, radiation, and relative humidity, among others. These non-precipitation factors are widely used to estimate potential evapotranspiration. However, this study focuses solely on pr. This is primarily because drought alleviation mainly depends on the effective management and allocation of available water resources. In future studies, the framework could be further improved by considering these non-precipitation factors. Another limitation concerns the range of pr increase scenarios considered in our experiments. As shown in Fig. S3, the long-term mean and standard deviation of monthly pr are approximately comparable, indicating that the maximum increase corresponds approximately to one standard deviation. Increases beyond 100 % would therefore represent stronger anomalies that are less frequently encountered and were not considered. Future studies could assess the drought alleviation response under such conditions. It is also essential to recognize that our study is a retrospective analysis based on historical MYDs. This study provides a framework for exploring attenuation efficiency within the early-stage response window of MYD development. However, determining whether a drought event would develop into an MYD at an early stage remains a major challenge (Xu et al., 2021; Gumus, 2023). This is due to the complexity of the factors influencing the development of droughts into MYDs. Furthermore, this study currently focuses solely on meteorological drought. Research on other types of drought, such as agricultural, hydrological, and socio-economic droughts, will also be a key focus in future studies. Future research should evaluate the applicability of our framework to other types of drought. To this end, it is also necessary to further optimize the framework by fully considering additional influencing factors (Tian et al., 2022; Wu et al., 2022b; Meresa et al., 2023).

6 Conclusion

This study applied a series of numerical experiments involving hypothetical precipitation gradient increases to investigate the attenuation efficiency for 351 historical MYDs during 1991–2020 across ten river basins in China, based on monthly PDSI calculated from daily meteorological data. Following these experiments, we developed the DAQM, a single-parameter model designed to quantify drought attenuation efficiency through the key coefficient, k, with t_optimal identified as the timing corresponding to k_max. The main conclusions can be summarized as follows:

(1) A nonlinear response relationship was identified between DS and hypothetical precipitation forcing applied incrementally along a gradient. Based on this relationship, the DAQM quantified attenuation efficiency and identified t_optimal, and showed strong agreement with both the mode-based and the small perturbation-based estimates, confirming its robustness and applicability for characterizing attenuation efficiency during MYD evolution.

(2) Our results revealed the existence of an early-stage response window during MYD evolution, corresponding primarily to t₁∼t₃. Within this window, attenuation efficiency is jointly influenced by the cumulative temporal advantage associated with earlier precipitation input and the higher numerical sensitivity of standardized drought indices under more severe drought conditions. Although t₁ remains the dominant t_optimal in most regions, accounting for 58.79 % of all grids, t₂ and t₃ are also non-negligible, accounting for 22.11 % and 11.06 %, respectively. These findings do not suggest that later precipitation input is generally more effective than earlier precipitation input. Instead, they indicate that within the early-stage response window, precipitation input during t₂ or t₃ may still exhibit attenuation efficiencies comparable to, or exceeding, that of t₁ under certain conditions. This implies that, during the early stage of MYD evolution, missing t₁ does not necessarily mean that the opportunity for relatively high attenuation efficiency has already been lost.

(3) DD is a key factor influencing attenuation efficiency within the early-stage response window. Longer DD was associated with a greater likelihood of t₂ or t₃ emerging as the t_optimal, evidenced by R² values of 0.526 and 0.578, respectively. When DD exceeds 29 months, the MYD attenuation efficiencies of t₂ and t₃ become comparable to that of t₁. Compared with the humid river basins in Southern China, the arid and semi-arid northern regions exhibited a higher probability that k at t₂ or t₃ would exceed that at t₁.

Overall, this study identifies and characterizes an early-stage response window during MYD evolution and further clarifies the competing mechanisms controlling attenuation efficiency within this window. With the increasing frequency and severity of MYDs in the future, our findings provide theoretical support and methodological guidance for improving the understanding of how precipitation input timing influences drought alleviation, thereby contributing to the development of more targeted and effective drought alleviation strategies across different regions.