We thank the reviewer for taking the time to carefully read our manuscript and provide constructive feedback. We value your efforts in helping us to make improvements to our paper. Below, we provide our responses to each comment and describe the corresponding revisions in the manuscript.

Our response is organized in the following order:

Reviewer Comment

Author Response

Paper Revision (if applicable)

Overall, I've found some confusion with the use of landslide "hazard" and "susceptibility" terms. I'd suggest checking the whole text and avoiding misuses of these terms. 

Because our modeling framework explicitly incorporates seismic and hydrologic forcing (e.g., peak ground acceleration and event-based recharge), the primary outputs of this study are landslide hazard maps. We have revised the manuscript to reflect this more clearly, and have removed or corrected instances where the terms were previously used interchangeably.

Was the Montecarlo approach used in the Landlab LandslideProbability component somehow constrained considering the ranges of parameters typical of the study area?

The Monte Carlo approach is constrained through parameter ranges that are informed by prior literature. 

To improve clarity, we have revised the manuscript (Line 210) to explicitly state that parameter bounds are defined relative to the mode and vary by parameter, following the Landlab LandslideProbability implementation, due to limited availability of region-specific reported minimum and maximum values. 

''Triangular distributions were used to represent parameter uncertainty because they provide a way to approximate variability around a central value when limited information is available. This approach is commonly used in landslide hazard modeling where site-specific measurements are limited but reasonable parameter ranges can be inferred from previous literature (e.g., Hammond et al., 1992; Selby, 1993; Strauch et al., 2018). Parameter bounds are defined relative to the mode and vary by parameter, following the Landlab LandslideProbability implementation.”

Was the landslide inventory somehow validated in the field? I'm asking this given that the event is relatively recent. 

While our initial landslide inventory relied on remote-sensing change detection methods (e.g., dBSI) over a large area, we recognize that reflectance-based classification can be affected by natural variability in soil, vegetation, and illumination, and thus provides limited direct validation of mapped landslides. Furthermore, this large-scale inferred inventory was not field-verified; only a subset (~20–30%) of the separately mapped coseismic landslides could be corroborated through field observations.

To improve robustness and provide a more defensible validation, we have shifted our approach to a terrain-based mapping method over smaller, high-priority areas within subcatchments that have pre- and post-event observations (Askerhan; Fig. 2) where landslides were concentrated following the February–March 2023 events.

In these areas, we performed detailed mapping of landslide source, transition, and depositional areas using high-resolution DEM differencing and slope analysis derived from lidar and/or photogrammetric data, allowing us to manually verify landslide locations and extents. This focused, hand-mapped approach enables a more defensible validation of our landslide susceptibility models, particularly for evaluating false negatives and model performance in capturing known initiation areas.

The revised, terrain-based landslide inventory does not alter the performance of the legacy model, indicating that its predictive skill is robust to improvements in mapping methodology. In contrast, the non-legacy model shows a notable improvement in performance, with AUC increasing to ~0.7, bringing it into alignment with the legacy model. As a result, both models now exhibit comparable discriminatory ability under the updated validation dataset. In addition to these quantitative results, we include qualitative, side-by-side comparisons of mapped probabilities, which show consistently higher predicted likelihoods in the legacy scenario. This pattern is physically consistent with the expected post-seismic weakening effect and supports the interpretation that inclusion of legacy processes improves hazard representation.

Regarding the sampling of non-landslide grid cells, why an equal number of landslide and non-landslide grid cells was selected?  I wonder if it would have been better to select a larger number of non-landslide cells than the same number of landslide inventory cells. This is a common approach in such types of analyses. 

Following the comment, we revised our sampling approach to better reflect the spatial distribution of mapped landslides. Rather than selecting an equal number of landslide (LS) and non-landslide cells (NLS), we now sample NLS cells in proportion to the total area of mapped landslide polygons. This ensures that the number of NLS samples is comparable to the total mapped landslide footprint, providing a more realistic assessment of model performance across the landscape.

This updated approach helps to align us with common practices in landslide susceptibility modeling (e.g., Strauch et al., 2018) and allows us to evaluate the model under conditions that better represent the imbalance between landslide and non-landslide areas.

How was the slope threshold of 15° used to filter out landslides (line 285) selected? 

In our previous analysis, we applied a slope threshold of approximately 15-20° to exclude transition and deposition zones and to minimize noise in the dBSI-derived landslide classification, consistent with physically based analyses indicating that slopes below this range are generally stable across the full spectrum of saturation conditions (Montgomery and Dietrich, 1994). We note, however, that in response to reviewer feedback, we have moved away from reliance on remote-sensing change detection and slope-based filtering, and instead now employ terrain-based mapping in a high landslide density area of our subcatchment, allowing for more accurate and defensible identification of landslide source areas. 

Line 385: Please add here details on the temporal resolution of the IMERG data used.

Line 385 has been updated to clarify the temporal resolution of the IMERG data. The revised text now reads: “Hydrological forcing is represented using daily-aggregated precipitation derived from the IMERG Version 7 product, which provides half-hourly global precipitation estimates."

Line 419: Please add here details on how the return period of the rainfall was calculated

Line 419 has been updated to clarify how the return period of the rainfall was calculated. The revised text now reads: "Within the study domain, max peak daily rainfall observed by IMERG during the AR event was 104.6 mm d-1, corresponding to an approximate 5-year return period estimated using Weibull exceedance probability applied to annual maxima of spatially maximum daily IMERG precipitation over the historical record (2000–2024; Fig. 8b)."

Some technical suggestions for the figures: 

Figs 1 and 2. An inset with the location of the study area would be useful

Insets showing the location of the study area have been added to Figures 1 and 2 to improve spatial context. These updated figures have been attached to this response.

Figs. 3-4. Please add details on the reference system (e.g. EPSG) of the maps 

Figure captions for Figures 3 and 4 have been updated to include details on the reference system. The revised captions now indicate that all data are projected in WGS 84 / UTM Zone 10N (EPSG:32610) for Landlab modeling. Example from Fig. 3:

"Illustration of model-data integration of remote sensing products into Landlab as fields: a) Copernicus GLO-30 DEM, b) USGS ShakeMap PGA, c) ESRI Sentinel-2 land use/land cover map, d) root cohesion derived from LULC in panel (c) using literature derived values (Strauch et al., 2018), e) MTA geologic map, and f) HiHydro Ksat. All data are projected in WGS 84 / UTM Zone 10N (EPSG:32610) for Landlab modeling."

Fig. 6. Please check the readability of the text in panels c and d. Moreover, please define all variables in the caption (nnl and nl are missing)

For Figure 6, the font size of the text in panels (c) and (d) has been increased and the lettering bolded to improve readability. Additionally, the caption has been updated to define all variables: “For each probability bin, the number of grid cells classified as legacy (n_l) and non-legacy (n_nl) is indicated in the figure.”

We also applied similar improvements to Figure 10, increasing font size and bolding letters in panels to enhance readability. The caption now clarifies the variables: “For each probability bin in (c), the number of grid cells classified as historic data based (n_h) and single-AR based (n_ar) is indicated in the figure."

These updated figures have been attached to this response.

We thank the reviewer for taking the time to carefully read our manuscript and provide constructive feedback. We value your efforts in helping us to make improvements to our paper. Below, we provide our responses to each comment and describe the corresponding revisions in the manuscript.

Our responses are organized in the following order:

Reviewer Comment

Author Response

Paper Revision (if applicable)

Discussion section must be separated and expanded. The interpretation of results, comparison with previous studies, and implications for hazard response are currently embedded in the Results section. A stand-alone Discussion is needed to: Explain the physical meaning of the observed reduction in critical slope angle under legacy conditions. Situate findings in the context of documented post-seismic landslide sequences (e.g., Wenchuan, Gorkha, Palu). Clarify how the proposed approach could inform rapid-response hazard assessment and early warning workflows. Explicitly discuss limitations, data uncertainties, and model assumptions.

We agree that a clearer separation and expansion of the Discussion improves the manuscript, and we have revised the structure accordingly by introducing a stand-alone Discussion section and separating it (line 709) from the Conclusions (e.g., 6. Discussion, 7. Conclusion).

We note that direct comparisons with previous studies were not extensively presented in the Results section; however, we have moved implications for hazard response into the Discussion, where they are more appropriately contextualized. We retained brief interpretation within the Results section to guide the reader, as these statements are concise and intended to provide immediate context for the presented outputs.

To address the specific points raised by the reviewer, we have expanded the Discussion as follows:

Physical interpretation of reduced critical slope angle: We added the following text to clarify the underlying mechanics (line 700): "This reduction in critical slope angle reflects the physical effect of seismic weakening on hillslopes, where shaking reduces slope shear strength through degradation of cohesion in soil and rock due to microfracturing, cracking, and regolith loosening. This lowers the factor of safety threshold for failure and allows previously stable slopes to fail at shallower angles (Brain et al., 2017; Keefer, 1984)."

Context within documented post-seismic landslide sequences and implications for hazard response: We expanded the Discussion to better situate our findings within prior work and to clarify the contribution of our approach. Specifically, we added (line 698): "Our framework leverages spatially distributed, gridded parameterization and a simplified representation of subsurface hydrology to efficiently estimate spatially distributed recharge and slope stability across large domains using minimal input data. This enables rapid, landscape-scale evaluation of evolving landslide hazard following an earthquake, providing a practical pathway for integrating earthquake legacy effects into operational hazard assessment and early warning workflows, particularly in data-limited regions."

In addition, we expanded citations to prior studies of post-seismic landslide activity to strengthen this comparison.

Model limitations, assumptions, and uncertainties: We have added a dedicated subsection (Section 6.1, Model assumptions and uncertainties) that explicitly outlines key limitations and sources of uncertainty in the modeling framework: "The purpose of this paper was to integrate coseismic and post-seismic landslide risk analysis which are often done separately in existing models. To serve this purpose we used parsimonious models like the infinite slope stability factor of safety model for landslide probability, steady-state subsurface flow for pore-pressure estimates, and a single-layer soil moisture accounting model for recharge. While the infinite slope approximation enables computationally efficient, landscape-scale probabilistic hazard estimation, it introduces uncertainty in areas with complex topography or heterogeneous materials, where failure surfaces may be non-planar, laterally constrained, or influenced by subsurface layering and preferential flow paths. As a result, modeled probabilities may underestimate hazard in deep-seated or compound landslides, or overestimate hazard on gently sloping, discontinuous terrain, emphasizing that our results reflect relative hazard and probabilistic susceptibility rather than precise predictions of individual landslides. Soil, root, and hydrologic parameters are parameterized locally using literature-derived mode values and Monte Carlo sampling to capture variability, including the use of triangular distributions to approximate uncertainty where site-specific data are limited. Important sources of uncertainty include heterogeneity in soil and vegetation properties, parameter estimation from gridded products, the scaling of hydraulic conductivity to account for lateral versus vertical anisotropy, simplifications in soil thickness parameterization, and the linear approximation of post-seismic weakening. Finally, steady-state hydrologic assumptions used in Landlab simplify computation but introduce uncertainty in local saturation estimates. These factors should be considered when interpreting model outputs and highlight opportunities for future model refinement."

We believe these revisions fully address the reviewer’s concerns and substantially improve the clarity, structure, and interpretability of the manuscript.

Clarify novelty and contributions. The manuscript should more clearly articulate what is new relative to prior compound hazard or post-seismic weakening studies. The innovation claims can be strengthened by emphasizing: (a) use of historical extreme rainfall records for pre-event hazard forecasting, and (b) quantitative mapping of legacy-induced slope threshold shifts.

We have updated the new Conclusion section to more clearly articulate the novel contributions of this study relative to prior work on compound landslide hazards and post-seismic landscape response.

"While prior studies have examined landslide triggering from individual drivers such as seismic shaking or rainfall (e.g., Strauch et al., 2018, Keefer, 1984; Iverson, 2000), and recent work has explored the concept of compound drivers (e.g., sequential earthquake and rainfall effects), several important gaps remain. This study makes two key contributions to advancing compound landslide hazard assessment. First, we integrate long‑term extreme rainfall climatology into a probabilistic landslide model to condition pre‑event soil wetness and antecedent susceptibility prior to seismic forcing. By using historical extreme rainfall records to account for antecedent soil wetness, our approach allows proactive mapping of areas likely to experience elevated landslide hazard, because these pre-conditioned soils align with the effects of future rainfall. Second, we introduce a quantitative mapping of legacy‑induced slope threshold shifts by relating post‑seismic strength reduction to peak ground acceleration within a Monte Carlo stability framework. This spatially distributed parameterization of legacy effects enables systematic comparisons of hazard modification due to landscape weakening following major earthquakes."

In addition, we have expanded the Discussion section to emphasize use of historical rainfall archives as hydrologic drivers of our landslide model:

“Finally, our results show that historical rainfall extremes, when combined with earthquake legacy effects, reliably reproduce the spatial patterns of rainfall-driven hazard predicted using the March 14-15, 2023 AR event data. The composite historical recharge model identifies areas of elevated soil saturation and susceptibility consistent with those highlighted in the single-event scenario. Grid-cell comparisons between the historical and single-AR P(F) maps reveal strong correlation, while slope distributions within probability bins differ by less than 0.4o on median critical slopes. These findings demonstrate that precomputed historical rainfall-driven hazard maps can effectively anticipate post-seismic landslide risk without requiring real-time climate data, providing a practical tool for proactive hazard assessment and planning.”

Model assumptions and uncertainties require further explanation.

We have added a new section (6.1 Model assumptions and uncertainties) at the end of the discussion to clarify key assumptions and associated sources of uncertainty in our modeling framework:

“The purpose of this paper was to integrate coseismic and post-seismic landslide risk analysis which are often done separately in existing models. To serve this purpose we used parsimonious models like the infinite slope stability factor of safety model for landslide probability, steady-state subsurface flow for pore-pressure estimates, and a single-layer soil moisture accounting model for recharge. While the infinite slope approximation enables computationally efficient, landscape-scale probabilistic hazard estimation, it introduces uncertainty in areas with complex topography or heterogeneous materials, where failure surfaces may be non-planar, laterally constrained, or influenced by subsurface layering and preferential flow paths. As a result, modeled probabilities may underestimate hazard in deep-seated or compound landslides, or overestimate hazard on gently sloping, discontinuous terrain, emphasizing that our results reflect relative hazard and probabilistic susceptibility rather than precise predictions of individual landslides. Soil, root, and hydrologic parameters are parameterized locally using literature-derived mode values and Monte Carlo sampling to capture variability, including the use of triangular distributions to approximate uncertainty where site-specific data are limited. Important sources of uncertainty include heterogeneity in soil and vegetation properties, parameter estimation from gridded products, the scaling of hydraulic conductivity to account for lateral versus vertical anisotropy, simplifications in soil thickness parameterization, and the linear approximation of post-seismic weakening. Finally, steady-state hydrologic assumptions used in Landlab simplify computation but introduce uncertainty in local saturation estimates. These factors should be considered when interpreting model outputs and highlight opportunities for future model refinement.”

The rationale for parameter distributions (triangular distributions for mechanical properties, scaling of hydraulic conductivity, soil thickness interpolation) should be discussed. A brief sensitivity analysis or statement about dominant sources of uncertainty would improve credibility.

We have expanded the manuscript to clarify the rationale for key parameter distributions and scaling assumptions, and to explicitly acknowledge their associated uncertainties.

Specifically, we added text in line 210 to justify the use of triangular distributions for mechanical hydrologic parameters: “Triangular distributions were used to represent parameter uncertainty because they provide a way to approximate variability around a central value when limited information is available. This approach is commonly used in landslide hazard modeling where site-specific measurements are limited but reasonable parameter ranges can be inferred from previous literature (e.g., Hammond et al., 1992; Selby, 1993; Strauch et al., 2018). Parameter bounds are defined relative to the mode and vary by parameter, following the Landlab LandslideProbability implementation.”

We further clarified the treatment of hydraulic conductivity (line 380), where scaling is applied to account for anisotropy in lateral versus vertical flow: “This scaling reflects the tendency for lateral hydraulic conductivity to exceed vertical conductivity in soils, which can enhance downslope drainage and reduce pore pressure buildup. A minimum value of 0.5 m d-1 was imposed to avoid unrealistically low conductivity values that would otherwise promote excessive saturation leading to slope instability. These adjustments introduce uncertainty where site-specific constraints on hydraulic anisotropy are unknown.”

In addition, we expanded the discussion of soil thickness parameterization (lines 375–376), including literature support for the adopted range and an explicit acknowledgment that this simplified representation does not capture fine-scale variability due to local controls such as lithology and vegetation: “This range is consistent with the reported soil depths on steep mountainous hillslopes and is intended to avoid unrealistically thin or thick soil thickness values that could bias model stability predictions (Strauch et al., 2018). We note that this parameterization introduces uncertainty, as soil thickness can vary substantially due to local controls (e.g., lithology, vegetation) that are not resolved in this formulation.”

Evaluation of the landslide inventory requires clarification. The rainfall-driven inventory derived from dBSI could be affected by vegetation phenology and bare-soil misclassification. The authors should add discussion of potential biases and, if possible, provide a few qualitative site-based validation examples.

While our initial landslide inventory relied on remote-sensing change detection methods (e.g., dBSI) over a large area, we recognize that reflectance-based classification can be affected by natural variability in soil, vegetation, and illumination, and thus provides limited direct validation of mapped landslides. Furthermore, this large-scale inferred inventory was not field-verified; only a subset (~20–30%) of the separately mapped coseismic landslides could be corroborated through field observations.

To improve robustness and provide a more defensible validation, we have shifted our approach to a terrain-based mapping method over smaller, high-priority areas within subcatchments that have pre- and post-event observations (Askerhan; Fig. 2) where landslides were concentrated following the February–March 2023 events.

In these areas, we performed detailed mapping of landslide source, transition, and depositional areas using high-resolution DEM differencing and slope analysis derived from lidar and/or photogrammetric data, allowing us to manually verify landslide locations and extents. This focused, hand-mapped approach enables a more defensible validation of our landslide susceptibility models, particularly for evaluating false negatives and model performance in capturing known initiation areas.

To assess the reliability of the dBSI-derived rainfall-triggered inventory, we compared it against our manually mapped terrain-based inventory in the Askerhan subcatchment. This comparison revealed that the dBSI inventory captured far fewer landslide points overall, and while some coincided with source areas, many identified pixels were located outside initiation zones, in transition or deposition areas. These discrepancies underscore the influence of vegetation phenology, bare-soil misclassification, and illumination effects on the dBSI method, and indicate that the remotely inferred inventory could not be reliably used for quantitative model validation. This evaluation directly motivated our shift to high-resolution, terrain-based mapping, providing a more defensible basis for validating landslide susceptibility predictions.

The revised, terrain-based landslide inventory does not alter the performance of the legacy model, indicating that its predictive skill is robust to improvements in mapping methodology. In contrast, the non-legacy model shows a notable improvement in performance, with AUC increasing to ~0.7, bringing it into alignment with the legacy model. As a result, both models now exhibit comparable discriminatory ability under the updated validation dataset. In addition to these quantitative results, we include qualitative, side-by-side comparisons of mapped probabilities, which show consistently higher predicted likelihoods in the legacy scenario. This pattern is physically consistent with the expected post-seismic weakening effect and supports the interpretation that inclusion of legacy processes improves hazard representation.

Sequence analysis requires stronger physical interpretation. The finding that “rainfall-before-earthquake” yields the highest hazard is interesting but currently under-discussed. The authors should relate this to known soil saturation effects on shear strength, crack propagation, and regolith weakening.

We have expanded the discussion to provide a stronger physical interpretation of the rainfall-before-earthquake scenario. This text has been incorporated into Section 6. (Discussion) following line 702 in the revised manuscript to provide a clearer mechanistic explanation for the high hazard observed in the rainfall-before-earthquake scenario.

“When examining the sequence of events, we found that the wet coseismic scenario produced the highest landslide hazard, reflecting the combined effects of antecedent soil moisture and seismic shaking. The enhanced landslide hazard observed in the wet coseismic scenario can be explained by established controls of soil moisture on hillslope stability. Increased antecedent rainfall raises pore water pressures, which reduces effective normal stress, in turn lowering the shear strength according to the infinite slope stability equation (Iverson, 2000; Montgomery and Dietrich, 1994). In addition, wet conditions may enhance crack propagation and reduce cohesion in partially weathered regolith, preconditioning slopes for failure under dynamic loading such as seismic shaking (Sidle and Ochiai, 2006; Iverson, 2000). Elevated saturation during the event also interacts with seismic shaking, allowing slopes that were marginally stable under post-seismic legacy conditions to fail when the critical factor of safety threshold is reduced by peak ground acceleration (Keefer, 1984).”

Some methodological details (e.g., component coupling, parameter tables) could be moved to Supplementary Material to improve readability.

We appreciate the reviewer’s suggestion. We have opted to retain the methodological details, including component coupling and parameter tables, in the main text to ensure transparency and reproducibility of the modeling framework.

Several figures could benefit from clearer contrast between high and low hazard classes.

In the revised manuscript, we have adjusted the color schemes of Figures 6, 7, and 10 to improve contrast between high and low hazard classes, enhancing the readability and interpretability of the hazard patterns. We have also improved readability for Figure 4. These updated figures have been attached to this response. 

The Conclusion should be more concise and emphasize implications rather than repeat results.

We agree with the reviewer that the Conclusion should be more concise and emphasize the broader implications of our findings rather than reiterate results. In the revised manuscript, we have streamlined the Conclusion to focus on novel contributions and important findings.

“While prior studies have examined landslide triggering from individual drivers such as seismic shaking or rainfall (e.g., Strauch et al., 2018, Keefer, 1984; Iverson, 2000), and recent work has explored the concept of compound drivers (e.g., sequential earthquake and rainfall effects), several important gaps remain. This study makes two key contributions to advancing compound landslide hazard assessment. First, we integrate long‑term extreme rainfall climatology into a probabilistic landslide model to condition pre‑event soil wetness and antecedent susceptibility prior to seismic forcing. By using historical extreme rainfall records to account for antecedent soil wetness, our approach allows proactive mapping of areas likely to experience elevated landslide hazard, because these pre-conditioned soils align with the effects of future rainfall. Second, we introduce a quantitative mapping of legacy‑induced slope threshold shifts by relating post‑seismic strength reduction to peak ground acceleration within a Monte Carlo stability framework. This spatially distributed parameterization of legacy effects enables systematic comparisons of hazard modification due to landscape weakening following major earthquakes. Additionally, our approach demonstrates the following points:

• Incorporating earthquake legacy effects improves landslide prediction accuracy at our validation sites for post-seismic rainfall-driven landslides, suggesting that existing rainfall-driven landslide hazard models can be updated for post-seismic periods with parsimonious parameterizations of earthquake legacy. 
• Across the full study domain, the dry coseismic model yields conservative performance that reflects dry antecedent conditions and provides a baseline for evaluating added effects of legacy and hydrologic forcing. Accuracy could likely be improved by calibrating soil saturation to February conditions typical to the area.
• Incorporating earthquake legacy effects leads to a ~13° reduction in median critical slope for landslide initiation compared to the non-legacy model, indicating that slopes considered stable before an earthquake may become highly susceptible afterward. This suggests that pre-earthquake hazard maps may significantly underestimate post-earthquake risk, reinforcing the need to update hazard assessments following major seismic events.
• The scenario where the AR event precedes the earthquakes shows the greatest landslide hazard, pointing to a potentially important hazard scenario not seen in the historical record and highlights the need for more research on how the timing of storms and earthquakes influences landslide risk.

It should be noted that validating hazard maps remains uncertain, particularly in remote mountainous regions where landslide inventories are often incomplete. An additional and more comprehensive evaluation of model performance would benefit from further field investigations to identify landslides or signs of instability that may have been missed in our satellite and mapping derived inventory.”