An integrated modeling method for rainfall-induced geohazard chains: coupling landslide mobilization and flow-state transition

1 Introduction

Climate warming is intensifying extreme rainfall in mountain regions, thereby increasing the likelihood of rainfall-induced landslides, debris flows, and flash floods (Westra et al., 2013; O’Gorman, 2015; IPCC, 2022). In steep catchments, intense rainfall can trigger not only shallow slope failures but also generate rapid hillslope runoff and channelized floods. Loose material released by landslides may be entrained, fluidized, and transported downstream by flowing water, allowing the hazard impact to extend far beyond the initial failure area. Events such as the 1999 Vargas disaster in Venezuela, Typhoon Morakot in southern Taiwan, and the 2023 extreme rainfall event in Beijing show that intense rainfall can generate clustered landslides, channel sediment supply, and downstream debris-flow or sediment-laden flood hazards (Larsen et al., 2002; Borga et al., 2014; Chen and Chuang, 2014; Ma et al., 2025). These cascading processes pose direct threats to mountain communities, transportation corridors, and downstream settlements and therefore require modeling tools that can support sustainable disaster-risk reduction and resilient land-use planning (Ortiz-Giraldo et al., 2023; Nong et al., 2026). Accordingly, the key challenge is not only to predict where slope failure occurs but also to understand how landslide material enters hillslope runoff and channelized flow and how it subsequently transforms into debris-flow-like motion or becomes diluted by flash-flood-dominated transport.

Modeling rainfall-induced landslide–debris-flow–flash-flood hazard chains requires defining the source conditions generated by slope failure. Physically based models have advanced from terrain-controlled stability analyses to grid-based transient infiltration models and coupled hydrological–geotechnical schemes, enabling more mechanistic prediction of shallow landslide initiation (Montgomery and Dietrich, 1994; Baum et al., 2008; Simoni et al., 2008). Probabilistic approaches have further been used to represent the influence of regional-scale hydrological and geotechnical parameter uncertainties on landslide prediction (Raia et al., 2014). However, identifying potentially unstable locations alone is insufficient for hazard-chain modeling. The spatial extent, thickness, volume, and timing of released material jointly determine the sediment supply and initial conditions for subsequent hillslope runoff, channelized flooding, and debris-flow transport (Fan et al., 2017; Tuganishuri et al., 2025). Source-area modeling therefore must link rainfall-induced local instability with a geometrically and volumetrically meaningful estimate of landslide material release.

After source release, the dominant uncertainty shifts from whether slope failure occurs to how the released material transforms and propagates. Landslide material may remain near the source area or, under sufficient water supply, fragmentation, and topographic confinement, evolve into a high-concentration debris flow connected to the channel network (Guthrie et al., 2010; Trujillo-Vela et al., 2022). During propagation, the moving mass may entrain bed or hillslope deposits, thereby modifying its volume, momentum, and resistance and ultimately affecting runout distance and impact intensity (McDougall and Hungr, 2005; Iverson and Ouyang, 2015; Kean et al., 2013). When channelized water becomes dominant, or the solid fraction decreases rapidly, high-concentration landslide–debris-flow material may also be diluted into a flash-flood-dominated, sediment-laden flow. Thus, landslide release, hillslope transport, debris-flow development, and channelized flood routing should be treated as interconnected components of a continuous hazard pathway rather than as independent modeling steps (van den Bout et al., 2022; Du et al., 2023; Tao et al., 2026). The central question addressed in this study is therefore how physically defined landslide-source conditions can be linked to the continuous transition among slide-like motion, debris-flow-like transport, and flash-flood-dominated sediment-laden flow.

This study develops a coupling-oriented framework for modeling rainfall-induced, landslide–debris-flow–flood disaster chains. The methodological contribution lies primarily in the coupling strategy and pathway diagnosis rather than in treating the hydrological, stability, CA, or routing components as isolated new submodels. Specifically, the framework converts hydrologically activated slope instability into landslide sources with explicit timing, geometry, and volume and then routes the released material as a single evolving solid–water mixture. This design enables continuous representation of slide-like motion, debris-flow-like transport, and flood-dominated sediment-laden flow, allowing rapid-dilution and sustained high-concentration pathways to be distinguished. Applied to the 2023 Beijing extreme rainfall event, the framework provides a physically interpretable basis for regional cascading-hazard assessment and preliminary hazard screening.

2 Rainfall-induced landslide source mobilization

This section defines how rainfall forcing is converted into a spatially explicit landslide source with release timing and volume. The formulation is organized around source mobilization rather than around separate hydrological and mechanical submodels. Rainfall first updates a compact set of subsurface states; these states modify pore pressure and suction; and source release is accepted only when local instability is hydrologically activated, geomorphically suitable, and patch coherent.

2.1 Hydrological states controlling slope activation

The hydrological component is not intended to resolve all catchment-scale runoff processes. It retains the transient subsurface states that most directly control shallow slope activation: wetting-front depth, perched saturated thickness, groundwater thickness, and effective saturation:where Hh is the hydrological state vector, Zf is wetting-front depth, Hp is perched saturated thickness, Hg is main groundwater thickness, and Se is effective saturation at grid cell x and time t. Accepted rainfall increases mobile storage above the wetting front, while recharge transfers water toward the lower saturated domain:

H h x , t = Z f , H p , H g , S e x , t , (1)

d M f / d t = I − R u , (2)

where M_*f* is mobile water storage (m), I is accepted infiltration (m·s⁻¹), and R_u is recharge to the lower soil domain (m·s⁻¹).

The wetting front advances according to the local moisture deficit so that antecedent wetness controls the timing and depth of hydrological activation:

d z f d t = I θ s − θ i z f , (3)

where θ_i and θ_s are antecedent and saturated volumetric water contents.

When mobile storage exceeds the remaining unsaturated capacity, a perched saturated layer is diagnosed; the main groundwater thickness is diagnosed from the groundwater-table position:

H p = M f − S u / n e , H g = H − z w , (4)

where n_e is effective porosity, H is local soil depth, and z_w is groundwater-table depth. To retain hillslope redistribution without introducing a fully distributed Richards equation model, lateral flow in the perched and main saturated domains is represented in compact Darcy–Dupuit form:

q p = − K s H p ∇ Φ p , q g = − K s H g ∇ Φ g , (5)

where qk is lateral flux, K_s is saturated hydraulic conductivity, H_k is the saturated thickness of domain k, and Φk is the hydraulic head. k = p denotes the perched saturated domain, and k = g denotes the main groundwater domain. This term updates Hh and redistributes locally generated wetting along hillslope gradients before evaluating stability. At the catchment scale, Equations 1–5 are evaluated at the digital elevation model (DEM)-cell level. The parameters appearing in these equations can be assigned from gridded or soil/lithology-class fields where available. Therefore, local soil depth, hydraulic conductivity, moisture deficit, and groundwater condition control cell-specific wetting-front advance, perched saturation, lateral redistribution, and pressure-head response before the stability calculation. In the notation used below, uppercase H refers to soil depth or saturated thickness quantities in the source module, whereas lowercase h refers to mobile mixture thickness in the transport module.

2.2 Hydrologically activated source criterion

Rainfall affects stability through two coupled pathways: it increases pore pressure at the candidate slip plane and reduces suction-related apparent cohesion. The equivalent pore pressure at the candidate slip plane is diagnosed from the hydrological state as Equation 6 where u_c is pore-water pressure at the candidate slip plane, ρ_w is water density, g is gravitational acceleration, β is slope angle, and H_c is the effective water-column thickness acting on the candidate slip plane. H_c is diagnosed from the main groundwater thickness H_g and from perched saturation H_p only when the perched layer reaches or hydraulically affects the candidate slip depth. Thus, the same hydrological state H_h controls both pore-pressure build-up and suction loss. The local factor of safety is then evaluated as Equation 7 where F_S is the dimensionless factor of safety; the terms C′, C_r, and C_ψ denote effective soil cohesion, root cohesion, and suction-related apparent cohesion (Pa), respectively; τ_r and τ_d are the resisting and driving shear stresses (Pa); z_c is the candidate slip depth (m); β is the local slope angle (°); γ_s is the soil unit weight (N m⁻³); and φ′ is the effective friction angle. A landslide source is not defined by a low-F_S cell alone. It must also satisfy hydrological activation, geomorphic source suitability, and patch-scale mechanical coherence as shown in Equation 8:where Ω_s is the source domain, F_c is the local stability threshold, A_h identifies cells in which wetting-front advance, perched saturation, or groundwater rise has become mechanically relevant, G_s denotes geomorphic source suitability, F_S,A is the patch-scale factor of safety, and F_A is the corresponding patch-scale threshold. This criterion converts rainfall-induced instability into a physically meaningful source domain rather than a scattered set of unstable pixels.

u c = ρ w g H c ⁡ cos 2 ⁡ β , (6)

F s = τ r τ d = C ′ + C r + C ψ + γ s z c ⁡ cos ⁡ 2 β − u c tan ⁡ φ ′ γ s z c ⁡ sin ⁡ β ⁡ cos ⁡ β , (7)

Ω s = x : F s x ≤ F c , A h x = 1 , G s x = 1 , F S , A x ≤ F A , (8)

2.3 Progressive source expansion and release volume

After initial source selection, progressive failure is represented by a slope-guided cellular automaton. Only cells on the active failure front can transfer destabilizing disturbance, and a neighboring cell is incorporated only when it remains hydrologically activated, near-critical, and sufficiently disturbed:

Ω f t + Δ t = Ω f t ∪ j ∈ N s ∂ Ω f t : A h , j = 1 ; F s , j ≤ F e n v ; p j t + Δ t ≥ p c , (9)

In Equation 9, N_s(∂Ω_f ^t) is the slope-guided neighborhood of the active failure front, A_h,j is the hydrological activation indicator, F_env defines the near-critical stability envelope, and p_c is the critical disturbance threshold. The volume supplied to the subsequent transport model is accumulated from the release depth over newly failed cells as shown in Equation 10:where Vol is released solid volume, d_r is release depth, and A is grid-cell area. The output of this module is therefore both a susceptibility map and a release schedule with explicit timing, source geometry, and solid volume.

Vol t = ∑ x ∈ Ω f t d r x , t Δ A , (10)

3 Released material routing and flow-state transition

Once landslide material is released, the methodological focus shifts from source activation to coupling the time-varying sediment supply with downstream flow-state transition. The released mass is represented as one evolving solid–water mixture. Landslide-like, debris-flow-like, and flood-dominated behaviors are not treated as independent layers or separate models but as states of the same mixture controlled by concentration, fluidization, maturity, resistance, and sediment exchange.

3.1 Solid–water mixture representation

The mixture is described by solid-equivalent thickness h_s, water-equivalent thickness h_w, total mobile thickness h = h_s + h_w, and depth-averaged velocity u. The solid fraction is expressed as Equation 11 where S _v is the solid fraction. In the diagnostic classification, a high S _v indicates slide-like motion, an intermediate S _v indicates debris-flow-like transport, and a low S _v indicates flood-dominated sediment-laden flow.

S v = h s h s + h w , (11)

The solid and water components share the same velocity, while the solid component exchanges mass with the bed and the water component responds to external water input or loss:

∂ h k / ∂ t + ∇ · h k u = S k , k ∈ s , w , (12)

S s = E s − D s , (13)

In Equations 12, 13, h_k is the equivalent thickness of component k, S_k is its source term, S_s is the solid source term, E_s is solid entrainment, and D_s is solid deposition. For k = w, S_w denotes net water input or loss. Concentration diagnoses the instantaneous mixture state, but it is not sufficient on its own because two mixtures with similar C_s values may have different internal structure, travel history, and deformation state.

3.2 Fluidization and debris-flow maturity

The solid fraction S_v describes the instantaneous composition of the mobile mixture, but it does not fully determine how the released landslide material is transformed internally. A dense mixture may remain grain-supported, become fluidized, or evolve into debris-flow-like motion depending on water enrichment, travel distance, and shear deformation (Iverson, 1997). To represent this path-dependent process, two internal state variables are introduced: the fluidization index κ and the debris-flow maturity index Ψ. The activation of each controlling process is described by a bounded function:

S X ; X o n ; X f u l l = min ⁡ [ 1 , ⁡ max ( 0 , X − X o n X f u l l − X o n ] , (14)

where X is a generic controlling variable, X_on is the onset threshold, and X_full is the full-activation threshold. This function allows each process to increase gradually from inactive to fully active.

The target fluidization state is calculated as Equation 15

κ * = 1 − 1 − F L 1 − F W 1 − F S h , (15)

where F_L, F_W, and F_Sh are bounded activation terms for travel-induced, water-induced, and shear-induced fluidization, respectively. They are evaluated using Equation 14 with travel distance L, water fraction W = h_w/h, and shear intensity represented by velocity or shear-power-related variables. This complementary form allows multiple mechanisms to jointly promote fluidization, while preventing any single factor from independently imposing a fully fluidized state.

The debris-flow maturity target is calculated only when the mixture remains sufficiently concentrated:

ψ * = G C G κ 1 − 1 − G T 1 − G S , S v ≥ S v , ⁡ min . (16)

For dilute flood-dominated mixtures, Ψ^∗ = 0 when S_v < S_v,min. In Equation 16, G_C represents the debris-flow concentration window, Gκ represents the requirement of sufficient fluidization, and G_T and G_S represent travel- and shear-supported maturity development, respectively. Therefore, κ describes the breakdown of the original slide-like structure, whereas Ψ identifies whether the fluidized material has developed into debris-flow-like transport. Both variables are advected with the mobile mixture, enabling continuous transition among slide-like, debris-flow-like, and flood-dominated states without prescribing separate moving layers.

3.3 Unified resistance and sediment exchange

The resistance term is written as a continuous blend of slide-like, debris-flow-like, and flood-dominated components:

a r = w L a L + w D a D + w F a F , (17)

In Equation 17, a _r is the effective resistance acceleration; a _L, a _D, and a _F are slide-like, debris-flow-like, and flood-dominated resistance components; and w_L, w_D, and w_F are concentration-dependent weights.

The three end-member resistance components are retained in compact form as where μ_slide and μ_flow are the friction coefficients for dense slide-like and fluidized states, respectively; τ_dil denotes the reference dilatancy or structural resistance of dense slide-like material (Iverson and George, 2014); τ_res represents the residual structural resistance retained under strong fluidization; and a _int accounts for unresolved internal shear and grain-collision losses in the depth-averaged formulation. The parameters μ_D and ξ_D vary between transitional and mature debris-flow values depending on the maturity index Ψ. The concentration-window factor GC confines the frictional transition to the debris-flow concentration range, whereas τ_reg, τ_vis, and τ_quad represent regularized yield, viscous, and shear-rate-dependent stress contributions, respectively.

a L = g ⁡ cos ⁡ θ μ L κ + 1 − κ τ d i l + τ r e s ρ m h + a int , (18a)

a D = g ⁡ cos ⁡ θ μ D + g u 2 ξ D h + τ r e g + τ v i s + τ q u a d ρ m h , (18b)

a F = g n 2 u 2 h 4 / 3 , (18c)

Sediment exchange is represented through mass-balance accounting terms rather than as a separate mobile phase.

E s = E L + E D + E F , D s = D c a p + D s e t , (19)

In Equation 19, E_L, E_D, and E_F denote slide-like plowing, debris-flow entrainment, and flood-driven hydraulic erosion, respectively; D_cap denotes capacity-controlled deposition; and D_set denotes low-mobility settling. This accounting separates the initial landslide release from secondary channel scour, reworked deposits, and permanent storage, which is essential for interpreting whether downstream high-concentration transport is sustained or rapidly diluted. The key parameters used in the coupled source–routing framework are summarized in Table 1.

TABLE 1. Key parameter settings for the coupled landslide source–routing framework.

4 Case study

4.1 Research area

The case study was conducted in a mountainous valley subset of Zhaitang Town, Mentougou District, western Beijing rather than across the entire administrative area. The area is characterized by steep relief, deeply incised valleys, and well-developed slope–channel connections, providing favorable geomorphic conditions for shallow landslide initiation and subsequent sediment delivery into channels (Figures 1a–e). The local DEM and post-event remote-sensing image show that shallow landslides were mainly distributed along steep hillslopes and valley-side slopes, where failed material can be rapidly connected to downslope drainage pathways.

The 2023 extreme rainfall event associated with Typhoon Doksuri was selected as the triggering scenario. To ensure temporal consistency between landslide-source generation and subsequent flow routing, an observed 95 h hourly rainfall sequence was used, with a cumulative rainfall of 469.0 mm and a maximum hourly intensity of 99.1 mm h⁻¹. This sequence includes antecedent rainfall, a low-rainfall interval, and the main rainfall period, thereby preserving the hydrological conditions required for time-dependent source activation.

The landslide inventory was refined from the publicly available shallow landslide inventory triggered by the July 2023 extreme precipitation event in Beijing. Post-event remote-sensing imagery was further examined to correct local omissions, refine landslide boundaries, and reinterpret landslide polygons within the local study area. This procedure produced an updated inventory of more than 400 shallow landslide polygons, which was used as the observational benchmark for evaluating the landslide-source simulation and subsequent slope–channel hazard-chain modeling.

4.2 Hydrological model benchmarking

Before applying the source-area model to the real rainfall event, a controlled hydrological benchmark was conducted on an idealized uniform hillslope to isolate the effect of hydrological model structure. The hillslope was 80 m long with a slope angle of 30°, a soil depth of 1.5 m, and a representative response depth of 0.9 m. A uniform rainfall event of 20 mm h⁻¹ lasting 18 h was imposed. The bucket model, Green–Ampt model, and the proposed dynamic front–perched layer formulation were compared using the same soil depth, hydraulic conductivity, water-content parameters, and wetting-front suction. This uniform benchmark was used only to isolate differences among hydrological model structures; the subsequent real-terrain application was implemented on the DEM grid, where terrain and source-related parameter fields vary spatially.

The three models produced different wetting and pressure-head responses under identical forcing (Figure 2). The bucket model showed the fastest storage response, with the slope-mean relative water-content state reaching approximately 0.84 at the end of rainfall, whereas the Green–Ampt and dynamic front–perched layer models reached approximately 0.67. In terms of pressure head, the bucket model produced a broad positive pressure zone near the response depth, with a maximum head of approximately 0.35 m, reflecting its column-integrated storage assumption. The Green–Ampt model mainly represented wetting-front advance and suction dissipation, with little positive pressure-head build-up. In contrast, the dynamic front–perched layer model generated localized positive pressure near the response depth while maintaining a more gradual water-content evolution. Mechanically, these differences imply different F_S responses: the bucket model tends to produce broad positive pressure-head development and thus a widespread reduction in effective normal stress; the Green–Ampt model mainly represents suction dissipation but produces limited positive pressure, leading to a weaker F_S reduction; the dynamic front–perched layer formulation generates localized positive pressure near the candidate slip depth, where suction loss and pore-pressure build-up can become mechanically effective for shallow landslide initiation.

The benchmark was further extended to the real DEM and the observed Beijing rainfall sequence, with the recession stage tracked to approximately 110 h. Two representative cells were selected from the interpreted landslide inventory and the non-landslide background area, respectively. All models remained nearly unchanged during the early low-rainfall period and responded rapidly after the main rainfall intensified at approximately 60 h (Figure 3a). The bucket model produced the strongest storage response, with the relative water-content state approaching 1.0 near the rainfall peak and then declining during recession, whereas the Green–Ampt model leveled off at approximately 0.84. The dynamic front–perched layer model produced an intermediate response and differentiated the two terrain locations, with the landslide cell reaching approximately 0.90 and the background cell peaking slightly lower and dissipating faster. This contrast was more pronounced in pressure head (Figure 3b): the Green–Ampt model mainly showed suction dissipation and remained slightly negative, the bucket model produced a rapid positive pressure increase up to approximately 0.57 m, and the dynamic model generated a higher peak at the landslide cell (∼0.44 m) than at the background cell (∼0.34 m), followed by post-rainfall recession. Because F_S is controlled by both pore pressure and suction-related apparent cohesion, this landslide–background contrast provides the hydrological basis for the subsequent CA source-mobilization simulation: cells with localized pressure-head amplification and sustained suction loss are more likely to enter the hydrologically activated, near-critical source domain.

4.3 Landslide-source simulation and evaluation

Using the locally refined landslide inventory as the observational benchmark, we evaluated the spatial performance of the cellular-automaton (CA)-based landslide-source model against Transient Rainfall Infiltration and Grid-Based Regional Slope-Stability Model (TRIGRS)-based rainfall-infiltration stability analysis. The interpreted inventory contains 1,575 landslide cells within the study area (Figure 4a). The CA model delineates 575 released-source cells that satisfy the coupled hydrological activation, instability, and source-selection criteria, representing compact initiation and release domains rather than the full spatial extent of interpreted landslide polygons (Figure 4b). These localized source cells subsequently develop into 2,133 final damaged cells through neighborhood interaction and topographic constraints (Figure 4c). The final damaged area captures the main landslide-prone belts while remaining more spatially selective than the stability-only screening result. By contrast, TRIGRS produces a broader minimum factor-of-safety field, with 4,363 cells classified as unstable under F_S < 1 and 4363 cells under the calibrated threshold of F_S < 1.05 (Figure 4d). These results indicate that the CA model represents a process-oriented transition from localized source initiation to spatially constrained damage expansion, whereas TRIGRS provides a threshold-based delineation of potential instability. This distinction motivates the subsequent analysis of CA-simulated release timing, sliding depth, and object-scale morphology.

To further characterize the process represented by the CA model, we analyzed the release schedule generated under the same 95-h rainfall sequence. The rainfall forcing has a cumulative depth of 469.0 mm, and the CA release timing is directly aligned with this observed sequence. The simulation produced 575 released cells, corresponding to a final released area of approximately 0.09 km², with release depth recorded at each time step. The temporal results show that landslide-source cells were generated episodically during the rainfall process and accumulated progressively into the final released area rather than being prescribed as a static failure mask (Figures 5a,b). The object-scale comparison further links this release process to landslide geometry and source-volume characteristics. The refined inventory contains 479 interpreted landslide polygons, while the CA simulation generated 398 released objects, of which 213 intersected the interpreted inventory within the matching buffer. The simulated objects are characterized by their area, elongation, and median release depth, providing the geometric and volumetric information required for subsequent debris-flow routing (Figure 5c). Thus, Figure 5 extends the spatial comparison in Figure 4 by showing that the CA model provides rainfall-aligned release timing, sliding depth, and object-scale morphology in addition to the final source-area pattern.

To assess the robustness of the CA-based source-area prediction under spatial parameter heterogeneity, we summarized 12 random-field scenario groups to test how source-related geotechnical and geomorphic heterogeneity, including soil thickness, strength-related parameters, spatial correlation structure, terrain amplification, and channel-related filtering, propagates through source selection and CA expansion rather than exhaustively sampling all hydraulic parameters. The stochastic analysis shows that the mean final-failure probability is spatially concentrated along the main interpreted landslide clusters, indicating that the predicted source areas are not randomly distributed over the hillslope domain (Figure 6a). At the object scale, 183 of the 479 interpreted landslide objects have high stochastic support (≥0.70), whereas 153 objects remain weakly supported (<0.10), revealing both the reproducible and missing portions of the inventory (Figure 6b). The process-stage metrics show that the F_S-only stage captures only a limited portion of the inventory, with an area ratio of 0.12 and an object-hit rate of 18.4%. After source selection and final expansion, recall increases to approximately 0.64 and the object-hit rate reaches approximately 65%, while precision remains approximately 0.52 (Figure 6c). Inter-group comparisons further indicate that the predicted area ratio is more sensitive to random-field settings than F1 score and object-hit rate, suggesting that the spatial correspondence with the interpreted inventory remains comparatively stable despite variations in predicted extent (Figure 6d). Overall, the stochastic results show that the CA model provides probabilistic support for a substantial subset of interpreted landslides, while its process-based filtering improves object-level coverage beyond local stability screening alone.

4.4 Landslide–debris-flow–flood transformation process

The previous section focused on identifying rainfall-induced shallow landslide-source areas. Once the released landslide material enters hillslope runoff and channelized flow, however, the process evolves from local slope failure into coupled water–sediment transport. This section therefore examines the short-term transformation from landslide release to debris-flow and flood-like transport following the largest release event.

The largest landslide release occurred at approximately 79.0 h, with an estimated released volume of approximately 7.44 × 10³ m³. Taking this event as the starting point, the system solid fraction exhibited a clear pattern of rapid increase followed by fast dilution within the subsequent hour. The system solid fraction rose briefly after the release, with the mean S_v reaching 0.16 at +5 min, but then decreased rapidly to 0.10 at +10 min and 0.02 by +60 min. The high-concentration tail showed a similar decay: the 90th percentile of S_v dropped from 0.65 at +5 min to 0.39 at +10 min and 0.06 at +60 min. This indicates that the release generated a short-lived high-concentration disturbance that was rapidly diluted during subsequent channelized transport.

Spatially, high-S_v cells were mainly distributed near the source area, adjacent footslope zones, and local convergent flow paths, whereas the main channel remained dominated by low-S_v flow. During the first 5–10 min, the high-concentration disturbance was redistributed along local drainage paths, representing the most evident stage of transition from landslide material to debris-flow-like transport. After 20 min, high-S_v areas contracted markedly, and the system gradually became dominated by channelized transport. By +60 min, high-S_v cells had largely dissipated, indicating substantial dilution of the initial release-induced disturbance (Figure 7).

To evaluate the model performance in representing outlet deposition and channel siltation, the remote-sensing interpretation, field photograph, and simulation results were compared (Figure 8). The interpreted results show that deposition was mainly distributed at the main channel outlet and along localized silted channel sections. The field photograph indicates that the typical outlet deposit thickness was approximately 1.2 m. In the simulation, the maximum mixture-flow depth, H_mix, was mainly concentrated along the main channel and outlet area, with local values reaching approximately 4 m, reflecting a short-lived high-intensity transport stage. In contrast, the final deposition depth was mainly concentrated at the outlet and local channel sections, mostly within 0–1 m, with an outlet depth of approximately 1 m. This is comparable to the field-observed thickness. Overall, the simulation reproduced the main spatial pattern and thickness magnitude of outlet deposition and channel siltation, supporting the interpretation of rapid dilution and localized deposition after short-term sediment transport.

After the landslide-derived material entered the channel system, the persistence of non-flood states was further evaluated through system-scale flow-state partitioning and event-scale dilution metrics. At the system scale, the mobile solid volume increased episodically in response to successive release events and peaked at approximately 88–91 h, approaching 3.0 × 10⁴ m³ (Figure 9a). Flood-dominated transport accounted for the largest fraction throughout the main stage. Debris-flow-like states increased transiently during periods of clustered releases, whereas slide-like states appeared only as short-lived peaks associated with individual release events, indicating that high-concentration states were difficult to sustain. Event-scale analysis further shows that non-flood responses decayed rapidly after the 19 major release events. For most events, T₅₀ ranged from several minutes to approximately 20 min, while some early events did not reach the T₂₀ threshold within the 60 min window, suggesting that both release magnitude and background channel flow modulated dilution persistence (Figure 9b). Event-aligned statistics show that the mean F_DF+SL decreased rapidly from 1.0 at release to approximately 0.5 at +10 min and further declined to approximately 0.3 after +20–30 min (Figure 9c). Overall, landslide releases generated short-lived, high-concentration pulses, but subsequent water-dominated transport rapidly diluted them, causing the system to evolve toward low-concentration sediment-laden flood flow.

After dilution of the high-concentration disturbance, bed erosion and deposition further regulated the subsequent solid supply. Cumulative erosion was mainly distributed along channels and local convergent flow paths, indicating that secondary sediment supply occurred primarily within concentrated transport pathways (Figure 10a). From the mass-balance perspective, the cumulative released solid volume during the main rainfall stage was approximately 3.29 × 10⁴ m³, whereas additional scour accounted for only approximately 7% of this volume. This ratio represents the simulated secondary-scour contribution under the local channel morphology and the capacity-limited sediment-exchange closure rather than a general entrainment ratio. Meanwhile, cumulative deposition reached approximately 2.13 × 10⁴ m³, suggesting rapid retention or deposition of much of the released material during transport (Figure 10b). The mechanism-based decomposition shows that additional scour was dominated by direct bed erosion, with limited contribution from reworked deposits (Figure 10c). The channel–hillslope partition further indicates that direct erosion increments were mainly concentrated in channel-like zones, while hillslope-like zones provided only limited supplements (Figure 10d). Therefore, landslide release controlled the initial solid supply, whereas subsequent channel scour acted mainly as a local secondary source, insufficient to sustain high-concentration transport.

5 Discussion

5.1 Source–water balance, threshold sensitivity, and pathway generalization

Recent studies increasingly emphasize that debris-flow persistence depends on both initial sediment supply and channel wetness, entrainment capacity, and the evolving distinction among streamflow, hyperconcentrated flow, and debris flow (Könz et al., 2024; Stammberger et al., 2024; Yang et al., 2024). In the present case, landslide releases generated short-lived high-concentration pulses, but these pulses were rapidly diluted by flood-dominated channelized transport. Event-scale metrics show that F_DF+SL decreased to approximately 0.5 within approximately 10 min and to approximately 0.3 after 20–30 min. Meanwhile, simulated additional scour contributed only approximately 7% of the released landslide volume, a case-specific and potentially conservative estimate reflecting localized channel erosion and the capacity-limited sediment-exchange closure. Channel erosion therefore acted mainly as a local secondary source rather than as a sustained sediment supply. These results indicate that the Beijing event mainly followed a rapid-dilution pathway, in which the released sediment pulse could not be maintained against dilution and deposition.

A closure-threshold sensitivity analysis was conducted to examine whether this pathway interpretation depends on empirical flow-state transition thresholds (Figures 11a,b). The perturbations affect the inferred partitioning among non-flood, debris-flow-like, and slide-like states, especially for the debris-flow-like and slide-like area–time metrics. However, the routing corridor and mobile-tendency response remain nearly unchanged, and the spatial Jaccard similarity remains high for most threshold-perturbed classifications. Therefore, the empirical thresholds mainly influence local flow-state classification and persistence rather than overturning the first-order rapid-dilution interpretation. The numerical threshold values are case-conditioned and should be recalibrated where sediment availability, channel confinement, rainfall forcing, or rheological properties differ. However, the pathway framework itself is transferable: systems with limited secondary sediment supply tend to evolve toward rapid dilution, whereas systems with stronger entrainment or sustained sediment replenishment may evolve toward sustained debris-flow or entrainment-maintained pathways in the S_v−I_fm state space (Figure 11c).

5.2 Pathway-dependent dynamic hazard and mitigation implications

The dynamic hazard maps provide a process-based view of how landslide-derived material modifies downstream impact intensity during the hazard-chain process (Figure 12). Unlike the source-area map, which only indicates where sediment is supplied, the dynamic hazard index reflects the combined effects of mobile mixture thickness h, velocity magnitude, and solid fraction S_v. It is therefore interpreted as relative hazard, rather than a risk, because exposure and vulnerability are not included. The results indicate that hazardous conditions were transient and pathway-dependent. The high-relative-hazard area increased from approximately 0.358 km² in the early post-release stage to 0.917 km² during the main transport stage, showing that landslide releases temporarily amplified downstream impact along connected flow corridors. However, the late-stage contraction of active and high-hazard zones suggests that this amplification was not sustained. This supports the rapid-dilution interpretation: the event was controlled by short-lived sediment pulses superimposed on flood-dominated transport rather than by a persistent debris-flow-dominated process. Such pathway-dependent hazard information can help identify slope–channel connection zones and downstream transport corridors that require priority monitoring or mitigation.

5.3 Model scope and limitations

The proposed framework is intended to resolve regional-scale hazard-chain pathways rather than local-scale internal flow physics. Its strength is that rainfall-induced release, channel routing, dilution, and sediment exchange are represented within one continuous modeling framework. From an implementation perspective, the DEM-cell source evaluation, slope-guided CA expansion, single-mixture HLLC/WB routing, and optional GPU/vectorized kernels support efficient scenario-based computation. The checkpoint design also allows release-time or rainfall-scenario branches to be updated from saved states, which is useful for wide-catchment screening and near-real-time early-warning interpretation. However, the transport module remains a depth-averaged effective description. It can distinguish rapid-dilution and sustained high-concentration pathways, but it does not explicitly resolve vertical phase separation, grain-size sorting, pore-pressure diffusion, dynamic hydro-mechanical effects, or evolving granular fabric (Liu et al., 2026).

This limitation is important because recent long-term debris-flow observations show that velocity, flow depth, discharge, bulk density, grain-size distribution, yield stress, and slurry viscosity all vary substantially among surges and directly control debris-flow mobility and deposition (Schöffl et al., 2025; Wei et al., 2025). Therefore, the simulated flow-state transition should be interpreted as a regional-scale estimate of hazard-chain persistence, not as a complete local rheological description of debris-flow interiors. Future constraints from flow depth, discharge, pore pressure, sediment concentration, grain size, channel-deposition observations, and exposed-element data would help refine the threshold between rapid dilution and sustained debris-flow transport and further extend the framework toward quantitative risk assessment when exposure and vulnerability data become available.

6 Conclusion

This study proposed a coupling-oriented modeling method for extreme-rainfall-induced geohazard chains by linking landslide-source mobilization, time-varying sediment supply, and downstream flow-state transition. The main conclusions are as follows.

• Hydrological benchmarking showed that different infiltration representations produced markedly different subsurface responses under the same rainfall forcing. Compared with bucket-type storage and Green–Ampt infiltration, the dynamic front–perched layer formulation better captured localized pressure-head development near landslide-prone locations, indicating that shallow landslide activation depends not only on total infiltrated water but also on whether pressure-head build-up becomes mechanically effective at the candidate slip depth.
• Applied to the 2023 Beijing extreme rainfall event, the CA-based source module generated rainfall-aligned shallow landslide release areas with explicit timing, depth, and volume information. Compared with stability-based screening alone, the process-based source selection and expansion improved object-scale correspondence with the interpreted landslide inventory, while the random-field analysis distinguished robustly supported source clusters from weakly constrained areas.
• The downstream simulation showed that released landslide material mainly formed short-lived high-concentration pulses rather than a sustained debris-flow-dominated process. During the main rainfall stage, the released landslide volume was approximately 3.29 × 10⁴ m³, whereas additional scour contributed only approximately 7% of this volume, indicating that channel erosion acted mainly as a local secondary sediment source.
• The event therefore followed a rapid-dilution pathway from landslide-derived sediment pulses toward sediment-laden flood transport. This suggests that rainfall-induced geohazard-chain evolution depends not only on whether slope material enters a channel but also on whether the released sediment pulse can persist against dilution, deposition, and limited entrainment. This process-based distinction can support preliminary cascading-hazard screening and sustainable disaster-risk reduction in mountain catchments exposed to extreme rainfall.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Junfeng Miao, University of Science and Technology Beijing, China

Reviewed by: Sartsin Phakdimek, Suranaree University of Technology, Thailand

Kuan Liu, Xi’an University of Science and Technology, China