Watershed-Scale Physics-Based Coupled Hydrological and Hillslope Stability Modeling of Rainfall-Induced Landslides: Back-Analysis and Forward Predictive Simulations
- Kassem, Mirna
- Advisor(s): Zekkos, Dimitrios
Abstract
Severe storms can trigger hundreds to thousands of slope failures in mountainous regions. Predictive slope stability models play a key role in improving community resilience to this widespread hazard. However, several limitations still constrain the prediction of rainfall-induced landslides at a regional scale. Existing rainfall-induced landslide predictive models vary along three general modeling aspects: (a) predictive approach (empirical or mechanistic); (b) size of the modeled area (single hillslope to global-scale applications); and (c) prediction type (probabilistic or deterministic). For mechanistic models specifically, two additional aspects are relevant: (d) the hydrologic model (steady-state or transient conditions); and (e) the landslide failure mode (one-, two-, or three-dimensional). While numerous predictive models have been developed, most struggle to balance mechanistic rigor, predictive accuracy, and computational efficiency for regional-scale assessments. Many rely on simplified mechanistic assumptions, while others use empirical approaches that lack any physical basis. Even when more rigorous mechanistic models exist, high computational costs, limited scalability, and sensitivity to input availability and quality affect their performance. This dissertation leverages advances in landslide and topographic mapping, physics-based modeling, and computational tools to enable scalable, robust rainfall-induced landslide hazard predictions.This dissertation introduces CRISIS (Coupled Regional Rainfall-Induced and Seismic Slope Instability Simulations), a new physics-based numerical modeling framework. CRISIS integrates a pseudo-three-dimensional slope stability approach with a hydrological model selected by the user, here implemented with the Parallel Integrated Hydrologic Model (ParFlow). This hydrological model simulates fully three-dimensional, transient groundwater flow alongside key surface processes, including overland flow, evapotranspiration, and vegetation dynamics. The CRISIS framework operates in two complementary modes: (1) back-analysis of mapped landslide inventories and (2) forward predictive modeling. Landslide inventories from past storm events are back-analyzed to estimate spatially variable shear strength parameters. These parameters are iteratively refined as additional landslide events are incorporated, enabling improved representation of subsurface variability across large regions where direct characterization is impractical. The resulting shear strengths can then be used in the forward modeling framework to predict the location, extent, depth, and timing of landslides triggered by future rainfall events. Leveraging parallel computing, CRISIS efficiently simulates large domains and is implemented as an open-source Python package to facilitate broad application. The model is applied to several watersheds in Puerto Rico impacted by Hurricane Maria (September 2017), which triggered more than 70,000 landslides and debris flows.A targeted sensitivity analysis was conducted to evaluate the influence of key model inputs, including hydrological formulations, initial groundwater conditions, hydraulic properties, and subsurface shear strength parameters, on predicted slope failures. Results demonstrate that representing groundwater flow in three dimensions is valuable for capturing realistic, spatially heterogeneous, and gradual failure patterns. In contrast, one-dimensional flow representations tend to produce abrupt or underestimated failures, particularly under moderate rainfall events. Spatial variability in initial groundwater tables results in a combination of top-down infiltration-induced failures and bottom-up saturation-induced failures, whereas assuming a spatially uniform initial water table leads to unrealistic synchronized slope responses. Antecedent rainfall further influences slope stability by reducing matric suction and eventually increasing positive pore pressures, which in turn accelerates the initiation of failure and enlarges the extent of instability. The relationship between hydraulic conductivity and rainfall intensity governs both the timing and spatial distribution of failures. From a mechanical perspective, cohesion exerts a dominant control on steep slopes, where small variations significantly alter landslide extent, while friction angle plays a secondary role. Additionally, high-resolution digital elevation models (DEMs) are critical for resolving fine-scale slope features, whereas coarser-resolution DEMs are sufficient for hydrological simulations, offering opportunities to reduce computational demand without compromising predictive accuracy. These findings highlight the importance of coupled hydrological-mechanical processes and guide key modeling decisions.Characterizing subsurface shear strength at regional scales remains challenging due to spatial heterogeneity associated with lithology, weathering processes, and limited subsurface observations. In this study, shear strength is interpreted not as an intrinsic material property, but as an effective, aggregate resistance that reflects the combined influence of soil and rock fabric, particle interlocking, degree of weathering, structural discontinuities, partial saturation, apparent cohesion from matric suction, minor cementation, root reinforcement, and transient pore-pressure conditions. To address this limitation, this dissertation leverages mapped landslide inventories to infer spatial patterns in shear strength through back-analysis. The CRISIS model was applied to a watershed in the Utuado municipality of central Puerto Rico, using a detailed inventory of 433 landslides triggered by Hurricane Maria. Landslide location, area, and volume were used to constrain these effective strength parameters. A key advancement is the incorporation of spatially and temporally varying pore water pressures simulated using ParFlow, enabling a realistic representation of hydrological conditions at the time of failure. This approach contrasts with conventional methods that often assume steady-state groundwater conditions or rely on simplified one-dimensional infiltration models. The back-analysis successfully reproduced observed landslides and yielded shear strength parameters consistent with Mohr-Coulomb representations for silty sand materials in the region. These values should be viewed as bulk expressions of slope-scale behavior, integrating frictional resistance, suction-dependent apparent cohesion, fabric effects, and weathering-induced weakening. Results reveal systematic variations in these effective parameters with depth and slope. These trends suggest a transition from friction-dominated to cohesion-influenced resistance at approximately 2.9 m depth, likely marking the boundary between regolith and less weathered bedrock. Hydrological conditions required to trigger failure also vary with geomorphic position, with higher pore pressures observed on gentle slopes and near valley bottoms and lower pressures on steep slopes and near ridgelines. These findings emphasize the coupled influence of topography, subsurface structure, and hydrology on landslide initiation.Beyond predicting landslide occurrence, CRISIS also provides insight into the mechanisms driving failure. While most regional models estimate aggregate metrics such as landslide area density or basic landslide characteristics, they do not typically distinguish between different failure mechanisms. In contrast, rainfall-induced landslides are classified into two primary mechanisms: (1) bottom-up, saturation-induced failure associated with rising groundwater tables, and (2) top-down, wetting-induced failure driven by infiltration processes, including downward-propagating wetting fronts and the formation of perched water tables above less permeable layers. Application of the model demonstrates that these mechanisms are strongly controlled by hillslope position and initial groundwater conditions. Bottom-up failures are most common on steep slopes near valley bottoms, where groundwater tables are initially shallow, whereas top-down failures dominate steep slopes near ridgelines, where groundwater tables are deeper. Failure depth is primarily controlled by subsurface cohesion, with lower cohesion values producing shallower landslides. Subsurface material properties, including both hydraulic and mechanical characteristics, significantly influence the frequency, timing, and type of failure. Additionally, slope failures due to pore pressure buildup at the colluvium-bedrock interface are associated with exfiltration of groundwater from the bedrock. This process requires (1) fractured, permeable bedrock underlying soil and saprolite, (2) spatial variability in soil thickness that allows bedrock outcrops to act as recharge zones, and (3) depth-dependent variations in hydraulic and mechanical properties. Although this mechanism has been observed in field studies, to our knowledge, it has not previously been modeled at the watershed scale. Measuring rainfall remains challenging in mountainous terrain and during extreme events, making it critical to understand how storm characteristics, such as intensity and duration, affect slope stability. Empirical predictive models are efficient and require minimal data but provide only binary outcomes, lacking information on the timing and extent of failures. In contrast, mechanistic models offer physics-based predictions but are computationally intensive and often depend on high-resolution data. To address this gap, we introduce a simplified, physics-based framework that is both computationally efficient and less dependent on extensive input data. This framework bridges the simplicity of empirical models with the robustness of mechanistic ones, enabling scalable and reliable Landslide Early Warning Systems (LEWS). It was developed using outputs from CRISIS simulations of 39 synthetic rainfall events, with varying total rainfall and duration. Although computationally intensive, the simulations are run only once per region and then used to build normalized predictive relations. The framework was applied to Hurricane Maria in the Utuado watershed. The predicted landslide area density closely matched both CRISIS predictions and field observations, with reduced input and computational cost. Additionally, the simulation outputs were used to derive a mechanistic rainfall Intensity–Duration (I-Dtotal) threshold, which aligned closely with a newly fitted empirical I-Dtotal threshold for Puerto Rico (1959-2024), further reinforcing the reliability of the model results.