Increased nonstationarity of stormflow threshold behaviors in a forested watershed due to abrupt earthquake disturbance

Extreme earthquake disturbances to the vegetation of local and regional landscapes could swiftly impair the former hydrologic function, significantly increasing the challenge of predicting threshold behaviors of rainfall–runoff processes as well as the hydrologic system's complexity over time. It is still unclear how alternating catchment hydrologic behaviors under an ongoing large earthquake disruption are mediated by long-term interactions between landslides and vegetation evolution. In a well-known watershed affected by the Wenchuan earthquake, the nonlinear hydrologic behavior is examined using two thresholds with intervening linear segments. A lower rising threshold (TH_r) value (210.48 mm) observed in post-earthquake local landslide regions exhibited a faster stormflow response rate than that in undisturbed forest and grassland–shrubland regions, easily triggering huge flash-flood disasters. Additionally, an integrated response metric pair (integrated watershed average generation threshold TH_g−IWA and rising threshold TH_r−IWA) with areas of disparate land use, ecology, and physiography was proposed and efficiently applied to identify emergent catchment hydrologic behaviors. The interannual variation in the two integrated hydrologic thresholds before and following the earthquake was assessed to detect the temporal nonstationarity in hydrologic extremes and nonlinear runoff response. The year 2011 was an important turning point along the hydrologic disturbance–recovery timescale following the earthquake, as post-earthquake landslide evolution reached a state of extreme heterogeneity in space. At that time, the TH_r−IWA value decreased by ∼ 9 mm compared with the pre-earthquake level. This is closely related to the fast expansion of landslides, leading to a larger extension of variable source area from the channel to neighboring hillslopes, and faster subsurface stormflow, contributing to flash floods. Finally, we present a conceptual model interpreting how the short- and long-term interactions between earthquake-induced landslides and vegetation affect flood hydrographs at event timescale that generated an increased nonstationary hydrologic behavior. This study expands our current knowledge of threshold-based hydrologic and nonstationary stormflow behaviors in response to abrupt earthquake disturbance for the prediction of future flood regimes.

1 Introduction

Understanding and measuring the hydrologic processes from local runoff generation mechanisms to larger watershed scales are difficult due to their complexity and nonlinearity (Farrick and Branfireun, 2014; Ross et al., 2021; Scaife et al., 2020). Previous researchers have made some efforts to identify the integrated, physically based hydrologic processes of the rainfall–runoff relationship in order to predict and simulate catchment runoff behavior under different conditions (Fu et al., 2013b; Gutierrez-Jurado et al., 2021). However, this cannot be generally applied due to uncertainties about water cycle processes, climate inputs, and complex physiographic boundary conditions. The observed hillslope- or catchment-scale threshold runoff response shows an emergent hydrologic pattern (Fu et al., 2013a; Ross et al., 2021; Wang et al., 2022; Zehe and Sivapalan, 2009) that could be used to identify key hydrologic signatures across different spatiotemporal scales (Ali et al., 2013). The hydrologic threshold is the critical point in time or space at which abrupt change in stormflow response occurs (Ali et al., 2013). Below the hydrologic threshold, a small amount of stormflow enters the adjacent channel, but significantly higher runoff magnitudes are generally observed above the threshold (Tromp-van Meerveld and McDonnell, 2006; Wei et al., 2020; Zehe et al., 2007). A unified threshold-based hydrologic theory that could potentially advance catchment hydrology was extensively discussed during the 2011 American Geophysical Union (AGU) Fall Meeting (Ali et al., 2013) and subsequently continuously developed (Ali et al., 2015; Ross, 2021; Ross et al., 2021; Scaife et al., 2020). Theoretical advancements in hydrology can support the development of appropriate algorithms for more efficient predictive hydrologic models.

The process of threshold behavior generally shows different nonlinear shapes of hydrologic response for a storage–discharge relationship (Ali et al., 2013; Wang et al., 2022), such as the hockey stick function, the step or Heaviside function, the Dirac function, and the sigmoid function. The transition from below-threshold to above-threshold behavior for different diagnostic shapes suggests several mechanisms of water retention and release in the watershed. In the literature, runoff behaviors with a hockey stick shape were found at the hillslope (Fu et al., 2013a; Tromp-van Meerveld and McDonnell, 2006; Wang et al., 2022) and watershed scales (Farrick and Branfireun, 2014; Scaife et al., 2020; Wei et al., 2020). For example, Farrick and Branfireun (2014) identified a threshold value of 289 mm for gross precipitation (P) and the antecedent soil water index (ASI) in a 3.15 km² forested catchment in Mexico. The threshold behavior presented two linear runoff responses, which were generally controlled by the subsurface stormflow mechanism. While the runoff behavior in the abovementioned catchment seemed to follow a hockey stick shape, above-threshold stormflow amounts showed high variability (Scaife and Band, 2017; Zhang et al., 2021a). This phenomenon possibly increased the uncertainty of prediction for higher stormflow amounts and flood disasters. Wei et al. (2020) proposed a rainfall–runoff relationship with multiple thresholds with intervening linear segments to reflect the initial streamflow activation and larger flood response. Understanding the change from slow to rapid stormflow response and to a larger flash-flood hydrograph is vital. However, a clear picture of the physical connotations of threshold behaviors associated with the generation and development of flash flooding is still missing (Wei et al., 2020).

Hydrologic threshold signatures at the catchment scale, as a new diagnostic tool, can effectively evaluate the long-term variation in stormflow response to forest recovery following natural disturbances (Ali et al., 2013; Scaife and Band, 2017; Wei et al., 2020). Natural disturbances in hydrology and their associated effects are summarized in Table 1, and they are also categorized into acute disturbances (ADs) and chronic disturbances (CDs). ADs, which are abrupt or sudden, such as earthquakes, wildfires, snow and ice, and volcanic activity (Table 1), tend to trigger the most drastic hydrologic response and alter the hydrologic regime after the disturbance. Comparatively, CDs (Table 1), which are more gradual and mostly affected by climate change (Hwang et al., 2018; John et al., 2022; Scaife and Band, 2017; Seidl et al., 2017), generally lead to a progressive reduction in the forest canopy without immediate destruction of the soil–root system and bedrock (Bladon et al., 2019; Hoek Van Dijke et al., 2022). Abrupt disturbance events could significantly alter the original landscape configuration and structure as well as the vegetation–soil system (Fig. 1), readily resulting in high peak flows and catastrophic flash-flood disasters (Arheimer and Lindström, 2019; Hoek Van Dijke et al., 2022). For instance, the famous Wenchuan earthquake on 12 May 2008 triggered nearly 2.0 × 10⁵ coseismic landslides (Fan et al., 2018; Xu et al., 2013), leading to the rapid and widespread destruction of the vegetation–soil structure and fragmenting rock mass structures (Cui et al., 2012; Zhang et al., 2021b). These disruptions can reduce the canopy interception and shallow soil water storage capacity, increasing the amount of throughfall precipitation reaching the soil surface and the stormflow magnitudes, thereby contributing to the flash-flood hydrograph (Zhang et al., 2018). After the abrupt disturbance, the exposed bedrock in the trailing edge of the landslides easily induces Hortonian overland flow, and the generated loose deposited material with high soil conductivities in the lower part of the landslides generally increases subsurface stormflow with the microporous flow (Mirus et al., 2017; Zhang et al., 2018). Such hydrologic behaviors are related to the quick runoff generation mechanism with a short lag time, resulting in higher runoff potential (Fig. 1). However, the hydrologic signature (e.g., soil water movement and stormflow generation) of risky landslides on steep hillslopes was not easy to capture. During the recovery processes, the number of earthquake-derived geohazards affected by large rainstorms led to unstable forest shrinkage and landslide expansions at long timescales in a forest-dominated mountainous watershed (Fig. 1). The unstable disturbances from endogenous (earthquake) and exogenous (rainstorms and concomitant hydrogeological hazards) origins remarkably increased the uncertainty in the assessment of the hydrologic regime from disturbance to recovery and flood risk management (Seidl et al., 2017). Previous studies have suffered from over-calibrated hydrologic models (Chiang et al., 2019; Maina and Siirila-Woodburn, 2019; Tunas et al., 2020) and a lack of understanding of runoff generation mechanisms in exploiting the effects of natural disturbance events on streamflow response. The efficient identification of nonlinear hydrologic behaviors in an earthquake-affected watershed is very significant to exploiting the formation mechanism of flash floods. Furthermore, understanding and clarification of the post-earthquake long-term dynamics of hydrologic threshold patterns are still urgently needed.

Figure 1 The long-term evolution of landslides in a disturbed watershed before and after the Wenchuan earthquake disturbance.

Table 1 Category and summary of natural disturbances in hydrology and the associated effects.

^a CD denotes chronic disturbance. ^b AD denotes acute disturbance.

Additionally, the scarcity of long-term hydrometeorological observation data, as well as the inaccessibility of post-earthquake roads, is a limitation to understanding the flood response driven by acute disturbance. In this study, the relationship between antecedent soil water storage, rainfall, and runoff at a 5 min interval in an earthquake-affected watershed was investigated. To understand the long-term variation in a hydrologic regime affected by an earthquake and their dominant controls, integrated hydrologic thresholds of precipitation and antecedent soil moisture with areas of disparate land use, ecology, and physiography were proposed at the watershed scale. The specific objectives of this study are as follows: (1) to examine the heterogeneity in hydrologic thresholds for undisturbed and disturbed lands, (2) to identify how subsurface stormflow and variable source area (VSA) affected by earthquake-induced landslides control the heterogeneity in the nonlinear physical processes from rainfall to runoff at the watershed scale, and (3) to use integrated threshold behaviors and the linear runoff response to gain insight into the long-term dynamics and nonstationarity of hydrologic behaviors affected by the interactions between post-earthquake landslides and vegetation evolution. Assessment of the streamflow regime and flood risk in a watershed affected by large physical disturbances could benefit from effectively identifying the hydrologic threshold signatures that mainly affect a watershed's streamflow response (Ali et al., 2013; Ross et al., 2021; Zhang et al., 2021a).

2 Data and methods

2.1 Study area

This study was carried out in the 78.3 km², forested Longxi River (LXR) experimental watershed, on the eastern Tibetan Plateau, China (Fig. 2). Forested land occupied 96.9 % of the watershed before the 2008 Wenchuan earthquake (Zhang et al., 2021b) and mainly consisted of open-canopy coniferous, dark coniferous, and broadleaf forests. After the earthquake, forested land decreased by 19.9 % (Zhang et al., 2021b). Thus, post-earthquake hydrogeological hazards, such as landslides and debris flows, could lead to an unstable landscape vegetation recovery trend in the area (Fig. 1), significantly influencing the stability of the hydrologic function and stormflow behaviors of the watershed, from rainfall to runoff (Zhang et al., 2021b).

Figure 2 The Longxi River (LXR) experimental watershed (78.3 km²) and its location on the eastern margin of the Tibetan Plateau (a, b). Panels (a) and (b) also show the area of the watershed affected by the 2008 Wenchuan earthquake and the detailed monitoring stations, mainly comprising rain gauges (R1–R9), gauging stations (S1–S6), and soil water probes (SW1–SW4).

The soil types mainly consist of Haplic Luvisols, Chromic Luvisols, Dystric Cambisols, and Haplic Alisols. Surface soil hydraulic conductivity is high for forested land, with values of 10–200 mm h⁻¹ (Zhang et al., 2021a). Subsurface stormflow generated at the soil–bedrock interface under heavy-rainfall conditions is a dominant runoff source contributing to flash flooding (Zhang et al., 2021a). The elevation in the region ranges from 870 to 3284 m a.s.l. (meters above sea level) with high relief. The strong orographic effect generally results in higher amounts of rainfall in steep mountainous watersheds (Fig. 2b), readily leading to large flood disasters. In the earthquake-affected regions, analyzing the disturbance–recovery processes of landscape vegetation will contribute to understanding potential long-term evolution in mechanisms of runoff generation and flash-flood disasters.

2.2 Hydrometric observations

Open-field precipitation from nine rain gauges was automatically recorded at a 5 min interval (Fig. 2a), and the flow discharge with a high flow velocity and water level during the flood hydrograph was monitored at a 5 min interval and calculated based on the hydraulic entropy method (Bahmanpouri et al., 2022; Chen, 2013; Moramarco et al., 2004; Zhang et al., 2021a). Volumetric soil moisture content (θ) was recorded at a 5 min interval using the soil water probes (Zhang et al., 2021a). Each probe, equipped with eight sensors at a 10 cm depth interval, was installed 80 cm below the surface in the soil profiles (Fig. 2a). The monitored probes SW1 and SW2 at upslope and downslope topographic positions were located in undisturbed forest and grassland–shrubland, respectively, whereas SW3 and SW4 were located at a disturbed landslide. The depth-equivalent antecedent soil water index (DASI) at the start of each rainfall event was obtained (Zhang et al., 2021a). The DASI is indicative of the initial shallow soil water storage (Wei et al., 2020) and is generally calculated from the eight-layer soil moisture measurements at each soil profile (Farrick and Branfireun, 2014):

$$\begin{array}{}\text{(1)}& \mathrm{DASI}=\sum _{i=\mathrm{1}}^n{\mathit{\theta }}_i\left(D_i-D_{i-\mathrm{1}}\right).\end{array}$$

Here, θ_i indicates the average soil content between the i and i−1 soil layer (cm³ cm⁻³), where i= 1, 2, 3, 4 … n and n indicates the number of soil layers below the surface for the monitored soil depth of 80 cm, and D_i indicates the soil depth at the ith layer (10, 20, 30, 40, 50, 60, 70, and 80 cm, D₀= 0). The index is utilized to exploit the effects of antecedent wetness on the magnitude of hydrologic thresholds and emergent behavior at the hillslope and watershed scales.

2.3 Definitions of storm events and hydrologic thresholds

Storms are defined as events with > 4 mm of precipitation that are separated by more than 6 h (Farrick and Branfireun, 2014; Penna et al., 2011). A total of 47 events were identified in this experimental watershed during the period from June to August of every year from 2018 to 2020, possibly filtering out the uncertainty in the assessment of hydrologic behaviors from seasonal variation in the forest canopy vegetation (Hwang et al., 2018). For each event, the stormflow (Q_q) is separated from the flood hydrograph using a proposed two-parameter recursive digital-filter method (Eckhardt, 2005; Zhang et al., 2021a). The catchment threshold behaviors between Q_q and the variable of the sum of the DASI + event precipitation amount (P) for each land-use type were quantitatively assessed using piecewise regression analysis (PRA), and the hydrologic threshold values of the sum of DASI + P and slope parameters from each linear segment of the PRA function were calculated (Oswald et al., 2011; Zhang et al., 2021a). The different break points and slope paraments of PRA might contribute to understanding the broad controls of shifts from slow to fast stormflow response and flash-flood disasters (Oswald et al., 2011; Scaife and Band, 2017; Zhang et al., 2021a). Uncertainty in the visual assessment of hydrologic thresholds is typically increased by nonlinear and complex stormflow behaviors (Detty and McGuire, 2010). However, the automatic identification of thresholds and linear slope parameters with a maximum likelihood approach (Muggeo, 2003) could be effective.

2.4 Determination of nonstationarity in pre- and post-earthquake threshold behaviors

Earthquake-induced landslides can destroy the soil–vegetation system, reducing the water storage of shallow soil and the canopy vegetation and leading to a change in the hydrologic thresholds of the sum of DASI + P for the disturbed landslide land-use type. The hydrologic threshold for a landslide is different from other undisturbed land-use types in the watershed. Meanwhile, as the destruction–recovery processes of post-earthquake landslides (Fig. 1), the mutual conversion of land-use types further influences the magnitude of water storage in the shallow soil and canopy vegetation for each land-use type, possibly altering the magnitude of the hydrologic threshold of the sum of DASI + P at the watershed scale. Herein, to clearly understand the long-term threshold evolution and the integral variation in emergent hydrologic behavior before and following the earthquake at the watershed scale, an integrated watershed average (IWA) index for the thresholds considering different land-use types was proposed to characterize the emergent watershed stormflow behaviors. The IWA index mainly considers the processes of runoff generation in the watershed's underlying surface based on the principle and framework of runoff potential for curve numbers (Deshmukh et al., 2013). The underlying surface mainly comprises the land-use types, the shallow water storage capacity and physical properties of soils at different locations, and bedrock types. These factors play vital roles in the processes of runoff generation. Another dominant water source for event precipitation in the atmosphere is also taken into account. Therefore, the index is mainly determined by the area contribution ratio of different land-use types (R_i), the shallow water storage capacity at different locations (DASI_i), and event precipitation:

$$\begin{array}{}\text{(2)}& {\mathrm{TH}}_{i-\mathrm{IWA}}\left(x\right)=\sum _{i=\mathrm{1}}^n{R}_i\times x_i=\sum _{i=\mathrm{1}}^n{\displaystyle \frac{a_i}{A}}\times \left({\mathrm{DASI}}_i+P\right),\end{array}$$

where DASI_i is the initial shallow water storage capacity of the ith land-use type in the underlying surface (mm), P is the event precipitation (mm), x_i is the runoff generation threshold or rising threshold for the ith land-use type (mm), a_i is the area of the ith land-use type (km²), A is the watershed area of the study region (km²), R_i is the ratio of a_i to A (%), and n is the number of land-use types. TH_i−IWA(x) is an integrated watershed average index for the thresholds at the watershed scale, with areas of disparate land use, ecology, and physiography. Waterbodies, buildings, and roads are assumed to be impermeable, and the initial storage capacity is assumed to be zero. Based on Eq. (2), TH_i−IWA(x) could be calculated to identify and compare the pre- and post-earthquake variation in threshold behaviors at the watershed scale, reflecting the integrated watershed hydrologic effects under the long-term interaction and development of the post-earthquake vegetation–hydrogeological hazards.

3 Results

3.1 Stormflow threshold behaviors

The collected event precipitation amounts (P) from a series of 40 large storm events (P > 10 mm) (Zhang et al., 2021a) showed great variation (from 16.4 to 263.9 mm). At larger P values, the stormflow amounts (Q_q) generally had very large variability, similar to the results reported by Farrick and Branfireun (2014). This phenomenon increases the uncertainty in the estimation of the rainfall–runoff relationship. No significant linear statistical relationship between DASI and Q_q was found at any of the monitored sites (p > 0.05, r²≤ 0.017; see our previous research, Zhang et al., 2021a). Once the DASI was combined with P, a more visually evident, nonlinear threshold behavior for the DASI +P and Q_q relationship (p < 0.001) was observed, with two hydrologic thresholds or break points for each location (Fig. 3), i.e., generation threshold (TH_g) and rising threshold (TH_r).

Figure 3 The piecewise regression analysis of event stormflow amount (Q_q) plotted against the sum of event precipitation amounts (P) and DASI for the (a) forest, (b) grassland–shrubland, and (c) landslide land-use types. Undisturbed forest and grassland–shrubland represent the pre-earthquake period, as reported by Zhang et al. (2021a), and the disturbed landslide land-use type represents the post-earthquake period. The red lines indicate the liner fit of the piecewise regression for the P+ DASI variable at the 95 % confidence level.

Additionally, significantly lower values of both TH_g (91.98 mm) and TH_r (210.48 mm) were observed for monitored landslide land (Fig. 3, Table 2). The threshold values of TH_g and TH_r for the landslide land-use type decreased by up to 28.84 and 43.86 mm, respectively, compared with average values in the undisturbed forest and grassland–shrubland land-use types. This indicates a lower post-earthquake threshold with response metric pairs that can lead to larger flash-flood disasters. The slope parameters from m_i1 to m_i3 were different by an order of magnitude (Table 2). A lower value of m_i1 mainly depicts a more gradual process of runoff generation. Observed larger storms in the third phase readily led to higher m_i3 values of > 1, indicating fast stormflow response during the flash-flood hydrograph. For different land-use types, lower m_i3 values occurred for monitored landslides. This was mainly owing to the deficiency in canopy vegetation and soil water storage capacity, highlighting the contribution of watershed storage during the first and second phases to the abrupt response of flash flooding in the third phase.

Table 2 Comparison of parameters used to assess the three linear threshold behaviors of the DASI +P and Q_q relationships at the 95 % confidence level.

^a m_ij indicates the values in the slope parameter of the PRA equations from the jth phase for the i land type (i may be forest, grassland–shrubland, or landslide and j= 1, 2, or 3, as shown in Fig. 3).
^b Collected data in the row were reported by Zhang et al. (2021a).
^c The correlation is significant at the 0.01 level (two-tailed).
^d SEE is the standard error of the estimate in multiple regressions.

3.2 Threshold pattern variation in an earthquake-affected watershed

Due to the fragmentation of the post-earthquake watershed's landscapes (Yunus et al., 2020; Zhang et al., 2021b) and the unevenness of monitoring locations for different land-use types (Farrick and Branfireun, 2014), the assessment of the hydrologic threshold behavior at the watershed scale is largely uncertain and nonstationary. The scarcity of hydrologic information before a large earthquake and the inaccessibility of disturbed post-earthquake regions by road increase the difficulty of timely monitoring in earthquake-affected watersheds. This limits our knowledge of how abrupt earthquakes affect the stormflow threshold behavior at the watershed scale.

The integrated watershed average (IWA) index for two thresholds was calculated using Eq. (2) to characterize the long-term changes in stormflow threshold behaviors before and following the earthquake disturbance. Significantly lower TH_g−IWA and TH_r−IWA values were found during post-earthquake periods (Fig. 4a), and the lowest values (109.34 and 250.72 mm, respectively) of both thresholds occurred in the coseismic phase. For 3 years (2008–2011) following the earthquake, the values of both thresholds remained low due to the unstable evolution of post-earthquake hydrogeological hazards (Shen et al., 2020; Zhang et al., 2021b). After the key turning point in the hydrologic disturbance–recovery process in 2011, following the earthquake, the thresholds recovered rapidly (within 3–5 years, 2011–2013, following the earthquake) and then gradually stabilized and approached their pre-earthquake levels. A significant negative correlation (p < 0.05) between the hydrologic thresholds and peak discharges was observed in this experimental watershed (Fig. 4b). Shortly after the earthquake, a significantly lower average TH_r−IWA value of ∼ 9 mm at the watershed scale increased the peak discharge and event flood volume by up to 22.58 % and 25.15 %, respectively (Table S1 in the Supplement). Figure 4a and b show the variation in the stormflow threshold behaviors and event flood response from 2007 to 2018. Four phases were observed during the hydrologic disturbance–recovery process, effectively and rationally predicting the flood regimes associated with stormflow threshold behaviors due to earthquake disturbance.

Figure 4 Changes in (a) observed stormflow threshold behaviors, including the integrated watershed generation threshold (TH_g−IWA) and rising threshold (TH_r−IWA), and (b) a large flood event response selected from Zhang et al. (2021b) during the period from 2007 to 2018 before and after the earthquake, including peak discharge (Q_p) and event flood volume (V, total flood flow in a single event). In phase 1 (2007–2008), TH_g−IWA and TH_r−IWA abruptly decrease and peak discharge rapidly increases. In phase 2 (2008–2011), TH_g−IWA and TH_r−IWA contain low values triggered by the overlapping of post-earthquake active geohazards. In phase 3 (2011–2013), TH_g−IWA and TH_r−IWA abruptly increase and peak discharge rapidly decreases. Finally, in phase 4 (2013–2018), the hydrologic variables gradually stabilize and approach the pre-earthquake level.

4 Discussion

4.1 Controls on threshold behaviors

Threshold behaviors and nonlinear stormflow responses at the watershed scale are useful to understand the generation and development of flash floods (Wei et al., 2020; Zhang et al., 2021a). They might also supplement the threshold-based hydrologic theoretical framework (Ali et al., 2013). To extend the application and efficacy of the derived TH_g and TH_r, three large flood events (19 August 2019, 15 August 2020, and 29 October 2020; see Fig. 5), which varied with respect to rainfall intensity at a 5 min interval (I_5min), event accumulative precipitation (EAP), and discharge (Q), were further considered. The TH_g value (120.8 mm) was observed during rapid streamflow change, and abrupt change in stormflow and flood response was readily stimulated at the threshold of TH_r using a value of 254.3 mm (Fig. 5). Such emergent behavior and signatures at the watershed scale demonstrated the presence of critical points in time and space where runoff response rapidly changed. The threshold index was efficiently verified via the application of the magnitude of two threshold values during the flood hydrograph (Fig. 5) and their concomitant hydrologic variation with the discharges before and following the Wenchuan earthquake (Fig. 4). However, we also acknowledge the limitation that only the dominant hydrologic process of runoff generation was considered, whereas the important hydrologic process of flow concentration was mostly ignored. In a future study, such metrics will be involved in the two hydrologic processes of runoff generation and flow concentration in order to more efficiently reflect the watershed's hydrologic behavior.

Figure 5 The generation and development of the flash-flood hydrograph, for the (a) 19 August 2019, (b) 15 August 2020, and (c) 29 August 2020 events, with the variation in 5 min rainfall intensity (I_5 min), event accumulative precipitation (EAP), and flow discharge (Q) based on the derived generation threshold (TH_g) and rising threshold (TH_r).

The nonlinear shape of storage–discharge relationships could reflect the runoff generation mechanisms underlying the retention–release processes of event water input at the watershed scale (Ali et al., 2013; Kirchner, 2009), efficiently diagnosing the inherent change in watershed nonlinear hydrologic behavior. Above the critical threshold value, a rapid discharge response could be observed in Fig. 3. The bedrock depression storage on the soil–bedrock interface (Fu et al., 2013b; McDonnell et al., 2021) and antecedent soil moisture storage (Cain et al., 2022; Fu et al., 2013a; Zhang et al., 2021a) are the main factors controlling the magnitude of the generation threshold (TH_g), influencing the initial emergent behavior of the rainfall–runoff process. Some common characteristics associated with process-based interpretation are highly permeable soils and low permeability of bedrock at the hillslope with steep slopes (Farrick and Branfireun, 2014; Ross et al., 2021; Scaife and Band, 2017; Tromp-van Meerveld and McDonnell, 2006). These properties generally lead to a significant soil–rock interface, readily triggering subsurface stormflow on the interface under heavy-rainfall conditions. The value of Q_q above the rising threshold was 2 orders of magnitude higher than that value below the generation threshold (Dickinson and Whiteley, 1970; Zhang et al., 2021a). At the TH_g threshold value, the bedrock depressions on the hillslope could be filled with water from rapid rainfall infiltration while water spilled over the undulating soil–bedrock interface (Fig. 6a, b) and generated higher streamflow (Fig. 6c). Once above the TH_r value, the VSAs close to the channel or impermeable surface under high-rainfall-intensity conditions showed a significant expansion (Fig. 6d). This will significantly increase the hydrologic connectivity of hillslope riparian–stream areas and readily facilitate the occurrence of catastrophic flash-flood disasters (Fig. 6e, f). The abrupt flow process was mainly affected by subsurface stormflow and Hortonian overland flow generation at the foot of a slope, and the fast runoff generation processes were subjected to a large storm size and higher-intensity precipitation rather than antecedent soil moisture.

Figure 6 Schematic diagram showing the changes in watershed (a, d, g) variable source area; (b, e, h) hillslope runoff process as cross sections; and flow discharge hydrographs at (c) the generation threshold (TH_g), (f) the rising threshold (TH_r), and (i) the TH_r affected by the earthquake-induced landslides, respectively.

Figure 6g shows the large expansion of VSA related to landslides, which generally destroy the crown canopy, litter layer, and soil–root system on hillslopes. These destructive processes alter the original landscape, resulting in the formation of bare rock, steep slopes, and the accumulation of loose materials (Chiang et al., 2019; Zhang et al., 2021b). For instance, the Wenchuan earthquake induced ∼ 2.0 × 10⁵ landslides that decreased the forest coverage by ∼ 30 % (Cui et al., 2012), thereby altering the infiltration–runoff processes and the contribution of Hortonian overland flow from the exposed bedrock at the trailing edge of the landslide as well as subsurface stormflow from landslide-generated loose deposited material to the flood hydrograph in the channel (Fig. 6h; Hong et al., 2010; Ran et al., 2015; Sidle et al., 2017). Zhang et al. (2021b) applied a hydrologic model and field observations to demonstrate that landslide-induced bare land can increase the runoff potential by > 10 % at the watershed scale compared with pre-earthquake conditions. Peak discharges increased by 22.58 %–367.42 %, and the time to peak was advanced by 25 min. The phenomena of increased peak discharge and advanced time to peak are almost consistent with the findings of Tunas et al. (2020) in the Bangga River basin, which was affected by the 2018 Sulawesi earthquake in Palu, Indonesia. The large physical disturbances triggered by earthquake-induced landslides can efficiently reflect the reduced water storage of shallow soil and the canopy vegetation. The post-disturbance stormflow thresholds of DASI +P can show a lower value (Table 2) and require less water input into the underlying surface to activate the runoff generation and flash flooding.

4.2 Nonstationary hydrologic behaviors triggered by an earthquake

Strong earthquake-induced landslides can severely fragment original landscape ecosystems; however, this occurs with spatially uneven distribution characteristics, such as the back-slope effect (Wang et al., 2019; Xu et al., 2011; Zhang et al., 2021b), the hanging-wall effect (Beyen, 2019; Huang and Li, 2009), and the distance effect along the earthquake fault zone (Xu et al., 2011). These spatially unbalanced conditions generally lead to apparent uncertainty and nonstationarity in the hydrologic response and runoff generation in earthquake-affected regions, severely affecting the recognition of the formation and development of flash floods. Zhang et al. (2021b) first illustrated the change in runoff response on both sides of the Longxi River watershed due to the back-slope effect. They found that the landslide area density of 13.76 % on the right bank during the coseismic period was approximately 3 times that (4.84 %) on the left bank, exhibiting an undulating but unpredictable disturbance–recovery process. This phenomenon highlighted the importance of the spatially uneven distribution and dynamic nonstationarity at timescales of earthquake-induced landslide patches for an accurate assessment of runoff generation and the dynamic evolution of catastrophic flash floods.

Interannual variation in rainfall–runoff relationships manifests a slow undulating increase in thresholds after the earthquake (Fig. 4a), predicting the short- and long-term changes and nonstationarity in stormflow behaviors at long timescales. These long-term, nonunique interannual thresholds indicated that stormflow behaviors at the watershed scale are closely related to catchment geophysical characteristics (Oswald et al., 2011) or climate (Graham and McDonnell, 2010), changing vegetation dynamics (Hwang et al., 2018), and natural disaster disturbances (Ebel and Mirus, 2014). Moreover, at event timescales, greater thresholds (TH_r) generally occurred during wetter growing seasons (Scaife and Band, 2017; Wei et al., 2020). Above the TH_r, catastrophic flash-flood disasters were readily triggered by summer high-intensity convective storms that possibly activated preferential soil flows, facilitating greater and faster stormflow generation in the shallow subsurface. This is associated with the lateral extension of hillslope–channel hydrologic connectivity (Farrick and Branfireun, 2014; Ross et al., 2021), abruptly increasing VSAs at the watershed scale.

Figure 7 Conceptual model explaining the long-term interactions between (b, d, f, h) disturbed hydrologic behaviors during the (a, c, e, g) formation–development–recovery process of the landslides after the Wenchuan earthquake.

By filtering out the influence of climate variables at event timescales, the effects of earthquake disturbance on flash-flood hydrographs were evaluated to determine how long-term interactions between earthquake-induced landslides and vegetation might alter interannual stormflow thresholds. Figure 7 presents a conceptual model of the long-term evolution and interaction between disturbed hydrologic behaviors affected by the unsteady formation–development–recovery processes of landslides after the earthquake. The Wenchuan earthquake largely destroyed the original landscape and vegetation–soil ecosystem, triggering hydrogeological hazards (Cui et al., 2012; Ran et al., 2015; Shafique, 2020; Zhang et al., 2021b). After the earthquake, larger landslides with exposed bedrock and loose deposited material can lead to the expansion of high-runoff-potential (high-RP) and VSA zones (Fig. 7a). This is related to Hortonian overland flow on exposed bedrock and subsurface stormflow, with microporous flow generated by loose deposited material from the landslides. These processes facilitate quick runoff generation, contributing to the flood hydrograph. The generation (TH_g) and rising (TH_r) thresholds associated with canopy vegetation interception and antecedent wetness conditions were significantly lower after the earthquake (Fig. 4a, Table 2), readily leading to a larger flood peak discharge and a shorter time to peak (Fig. 7b). Due to the complex geohazard processes, the destruction of the vegetation–soil ecosystem and the widened channel further expanded the high-RP and VSA zones and enhanced the hydrologic connectivity structure shown in Fig. 7c (Moreno-de-las-Heras et al., 2020). The value of TH_r related to large flash floods drastically decreased (Fig. 7d), generating higher peak discharge and less lag time. After the key turning point in 2011, the rapid recovery of landslide areas and the vegetation–soil ecosystem was observed from 2011 to 2013 (Fig. 7e, f). During this period, the RP and VSA zones rapidly expanded and the hydrologic threshold behaviors quickly recovered and improved. Our findings emphasize the importance of earthquake-induced landslides and vegetation dynamics with respect to the event-scale and long-term stormflow response, but more research is required to fully understand the disturbed hydrologic behaviors and threshold-based hydrologic theoretical framework.

4.3 Challenges in identifying nonunique threshold behaviors disturbed by large disaster events

Abrupt disturbance events generally disrupt the hydrologic storage and functional connectivity at both spatial and temporal scales (Ebel and Mirus, 2014), readily leading to nonunique thresholds in the nonlinear behaviors of rainfall–runoff processes (Arheimer and Lindström, 2019; Scaife and Band, 2017). Uncertainties and challenges in characterizing the nonunique threshold behaviors from disturbance hydrology perspectives still need to be clarified.

The lack of detailed observation data about hydrologic fluxes and subsurface storage dynamics before and after an abrupt event is a significant limitation to understanding the variation in flow pathways and runoff generation mechanisms (Chiang et al., 2019; Ebel, 2020; Farrick and Branfireun, 2014). Thus, it will be difficult to accurately identify the presence and the form of hydrologic threshold signatures that influence a watershed's streamflow response, possibly leading to assessment uncertainties with respect to the flood hydrologic regime.

The watershed's spatially broken patchiness and dispersion brought on by the occurrence of sudden events are very evident (Sidle et al., 2017; Wang et al., 2019; Zhang et al., 2021b), but the methods to quantitatively assess the variable functional hydrologic connectivity associated with runoff generation mechanisms are still lacking (Beiter et al., 2020; Bracken et al., 2013). This limits our understanding of the effects of the post-earthquake hydrologic responses and long-term variable threshold behaviors on flooding in the disturbed regions.

Heterogeneity in the regional plant community composition and subsurface critical zone thickness related to the water storage capacity (García-Gamero et al., 2021; Hahm et al., 2019; Shangguan et al., 2017) generally leads to different runoff generation mechanisms and nonstationary threshold behaviors in different climate zones. A reasonably generalized ecohydrologic zoning with some organizing principles or frameworks could be a better road towards a unified threshold-based hydrology theory, as proposed by Ali et al. (2013), further facilitating cross-site syntheses and validation.

5 Conclusions

This study offered insight into the complex interactions between hydrologic processes in a small, disturbed, forested experimental watershed by revealing relatively simple, generalized, nonlinear runoff behaviors. An integrated response metric pair was applied to identify stormflow threshold behaviors and evaluate the long-term threshold dynamics following the Wenchuan earthquake disturbance. Lastly, we revealed that the subsurface stormflow and VSA are dominant controls on the dynamics of threshold behaviors before and following an acute disturbance, rather than a chronic disturbance. In summary, the key findings of this work are as follows:

• (1)

  Lower values for both the generation and rising thresholds derived from nonlinear stormflow behaviors generally occur in disturbed landslide regions, more easily leading to large flash-flood disasters.

• (2)

  The dynamics of the two thresholds, via a novel integrated watershed average index, can characterize the hydrologic disturbance–recovery process before and after an abrupt earthquake.

• (3)

  The runoff generation mechanisms of subsurface stormflow and VSA mainly control the nonstationarity in threshold behaviors and linear stormflow response. This is related to the largely spatially heterogeneous distribution of abrupt earthquake-induced landslides and the temporally undulating recovery of disrupted landscapes and vegetation–soil ecosystems. The abrupt hydrologic disturbance differs from the nonstationarity in vegetation–climate interactions, leading to chronic variable thresholds.

This study emphasizes the importance of integrated stormflow thresholds as a diagnostic tool to effectively characterize abrupt variation in emergent catchment patterns and the broad shift from slow to fast flood response, particularly with temporal nonstationarity in the long-term interactions between abrupt earthquake-induced landslides and vegetation evolution. This work contributes to the mitigation and adaptive strategies for unpredictable hydrologic regimes and flash-flood disasters triggered by abrupt natural disturbances.