The dominant role of latent spatial structure in landslide susceptibility

Aristizábal, Edier; Lombardo, Luigi; Korup, Oliver

doi:10.5194/egusphere-2026-1749

Preprints

https://doi.org/10.5194/egusphere-2026-1749

Preprints

18 May 2026

| 18 May 2026

The dominant role of latent spatial structure in landslide susceptibility

Edier Aristizábal, Luigi Lombardo, and Oliver Korup

Abstract. Landslide susceptibility mapping is essential for risk management and mitigation. Traditional multivariate models have achieved high nominal accuracy in predicting landslide occurrence, but the underlying methods rarely account for spatial dependence in the input data. We test how spatial Hierarchical Generalized Linear Models (HGLMs) and spatial autoregressive models might enhance both accuracy and reliability of landslide susceptibility. We estimate the frequency of landslides in a catchment, drawing on a catalogue of 10,837 landslides and several predictors (e.g., rainfall, elevation, slope) in the northern Colombian Andes. Our HGLMs integrate the effects of spatial dependency and heterogeneity through Markov Random Field (MRF) models, i.e. Intrinsic Conditional Autoregressive (ICAR), Besag-York-Mollié (BYM), and Leroux models. Our results show that landslide frequency is significantly influenced by spatially-dependent unobserved factors, likely representing contiguous geological formations or shared soil properties, which we therefore incorporate as latent variables. Even after accounting for known covariates, the HGLMs capture spatial dependencies that non-spatial models fail to address. Incorporating spatial structure in the data improves model performance, judging from model selection metrics such as the Deviance Information Criterion (DIC) or the Watanabe–Akaike Information Criterion (WAIC). By accounting for latent spatial effects, spatial HGLMs produce smoother and more reliable susceptibility maps. This approach overcomes a key limitation of traditional models: the underestimation of landslide frequency in high-density areas where unobserved, spatially-structured factors are most influential. Our findings highlight the importance of integrating spatial dependence and heterogeneity in landslide susceptibility models to achieve enhanced predictive performance and reliability.

Received: 27 Mar 2026 – Discussion started: 18 May 2026

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.

Download & links

Edier Aristizábal, Luigi Lombardo, and Oliver Korup

Status: final response (author comments only)

RC1:
'Comment on egusphere-2026-1749', Anonymous Referee #1, 27 May 2026
General assessment
This manuscript addresses an important methodological question, namely the role of spatial dependence in landslide susceptibility modeling and the use of hierarchical spatial models. The comparison of non-spatial and spatial models is potentially relevant for NHESS readers, and the manuscript is clearly structured and generally well written. However, I cannot recommend publication in its current form. My concerns are conceptual as well as technical: (1) The study operates at the level of terrain units that is not appropriate for landslide susceptibility modeling, which undermines model interpretation and relevance of findings. (2) The study conflates statistical concepts, potentially over-interpreting and over-stating its findings. I therefore recommend rejection.
Major concerns
1. Terrain mapping units are inadequate for susceptibility modeling
The manuscript models landslide susceptibility using hydrographic catchments of approximately one hundrd to a few hundred km² as terrain mapping units. This support is inadequate for a process known to be governed by slope-scale conditions. Yet the predictors are aggregated to catchment means (mean slope, mean elevation, rainfall summaries), suppressing the variability that controls slope failure. This severe change-of-support problem likely dilutes predictor-response relationships by construction, resulting in inflated residual variation. In my view, this issue fundamentally undermines the manuscript’s central conclusion.
2. Over-interpretation of latent spatial effects
The manuscript repeatedly interprets latent spatial effects as evidence for contiguous lithology, soil properties, tectonic domains, or shared physical controls. However, the CAR formulation only encodes adjacency among catchments. The model demonstrates that neighboring polygons have similar residual structure; it does not identify geology or tectonic controls. No attempt is made to relate latent spatial structure to known geological or other structures in this vast mountain area, and no attempt is made to use such (even coarse) information to construct additional process-related covariates. The physical interpretation is therefore speculative and substantially overstated.
The use of three large hydrographic basins as higher-level units is particularly weakly justified. No evidence is presented that these basins correspond to meaningful geophysical groupings for slope stability, and there is no general reason to expect “tectonic domains” to align with hydrographic catchments.
3. Model interpretation and diagnostics are insufficient
The manuscript repeatedly makes strong claims (“reliable”, “accurate”, “better estimates”, “overconfident predictions”) without clearly distinguishing precision and bias of coefficient estimtors, predictive performance, and goodness of fit. Residual diagnostics also raise concerns. Pearson residuals appear extremely large (up to 20-30), suggesting major lack of fit or overdispersion (perhaps related to points 1 and 2 above), yet this is not discussed.
4. Predictor selection and process representation are insufficiently justified
The predictor set is extremely limited (mean slope, mean elevation, rainfall summary). Relevant predictors commonly used in susceptibility modeling (e.g., land cover/land use, lithology, anthropogenic disturbance such as roads, topographic wetness index, and of course local slope angle) are absent without justification. Predictor selection (correlation analysis and PCA) is insufficiently described, and the manuscript does not justify linear predictor-response assumptions despite well-known nonlinearities, particularly for slope angle.

Additional comments
Critical information on the landslide inventory and the predictor variables is missing, e.g. concerning completeness, consistency in digitizing landslides, homogeneous availability of multitemporal satellite imagery. Size distribution information would help to assess completeness for small landslides.

Clarify coefficient interpretation across models. M1 appears to use unstandardized predictors and rate ratios, whereas M2–M5 report posterior means after standardization. Table 1 therefore appears internally inconsistent and coefficients are not directly comparable.

Report confidence intervals and credible intervals consistently and prefer interpretable effect sizes, e.g. 100 m elevation difference not 1 m.

Moran’s I should be defined in Methods, not introduced implicitly in Results/Figure captions.

The Results section contains methodological material (e.g., model definitions) that belongs in Methods.

The Introduction overstates hierarchical spatial autoregressive models as the primary solution while overlooking alternative approaches (e.g., spatial smoothers, autologistic models, spatial bootstrap approaches), some of which have already been used in the landslide literature.

The Discussion does not discuss the findings in the context of the literature.

Additional minor comments
L11 "model performance" - or rather "goodness of fit" in the usual statistical terminology; "performance" suggests estimation on independent test set

L22-23 "cross-sectional" and "panel data" are terms that are not commonly used in the environmental statistics literature; prefer more descriptive alternatives? e.g. longitudinal or simply temporal or time series?

L24 TMU is not a concept from spatial statistics

L39 "model residuals" - definition not so obvious in the context of classification problems

L41 "interactions" - concept not obvious in this context; unobserved confounders are not interactions - not statistically and not otherwise

L44-51 "statistical oversight" is very strong wording, considering that most susceptibility analyses do not interpret model coefficients in detail and do not perform statistical inference.
Citation: https://doi.org/10.5194/egusphere-2026-1749-RC1
- AC1: 'Reply on RC1', Edier Vicente Aristizábal Giraldo, 06 Jun 2026
  
  Major Concern 1: Terrain mapping units are inadequate for susceptibility modeling
  We respectfully disagree with this assessment. While landslides are inherently slope-scale processes, this does not preclude their analysis at coarser spatial units. The choice of analysis scale is governed by the study objectives, data availability, and the desired spatial resolution of the output map. Regional and national-scale susceptibility assessments, which necessarily operate at resolutions far coarser than the individual slope, are well established in the literature (Guzzetti et al., 2005, Geomorphology, 72:272–299; Stanley & Kirschbaum, 2017, Natural Hazards, 87:145–164; Nowicki Jessee et al., 2018, JGR: Earth Surface, 123:1351–1369). At those scales, predictors such as geology, mean slope, and climatic variables effectively serve as regional proxies for local susceptibility conditions.
  Critically, our dependent variable is not the binary presence/absence of a landslide on a slope, which would indeed require slope-unit or pixel resolution, but rather the count (and implicitly density) of landslides per catchment. This Poisson count formulation is specifically designed for aggregated data and is theoretically consistent with catchment-level predictors (Lombardo et al., 2020, Earth-Science Reviews, 209:103318). Mean catchment slope and elevation are appropriate summaries of terrain conditions that drive bulk susceptibility at this scale, a design choice analogous to using basin-wide lithology fractions or mean annual rainfall in regional hazard assessments. We will add a paragraph in the Methods section explicitly justifying the catchment unit in the context of the Poisson count model.
  
  Major Concern 2: Over-interpretation of latent spatial effects
  We partially agree and will revise the Discussion accordingly. The CAR/ICAR/BYM formulations encode spatial adjacency structure in the residuals; they do not identify specific geological or tectonic domains. Our Discussion text listing "geological characteristics, tectonic settings, and anthropogenic effects" was intended to enumerate plausible sources of unobserved spatial structure, not to assert that the model identifies them. We will rephrase this to make the distinction explicit, in line with standard practice in spatial hierarchical modelling (Besag et al., 1991; Lawson, 2018, Bayesian Disease Mapping; Morris et al., 2019, Spatial and Spatio-temporal Epidemiology, 31:100301).
  Regarding the three basins: the Atrato, Cauca, and Magdalena drainage systems correspond to distinct structural and lithological domains of the northern Colombian Andes, the Pacific-facing Western Cordillera (Atrato), the inter-Andean Cauca depression flanked by volcanic and metamorphic lithologies, and the Magdalena valley between the Central and Eastern Cordilleras, making them physically motivated grouping units for spatial heterogeneity at the upper level. We will substantiate this choice with supporting references in the revised manuscript.
  
  Major Concern 3: Model interpretation and diagnostics are insufficient
  We agree and will revise accordingly on both points.
  
  Major Concern 4: Predictor selection and process representation are insufficiently justified
  We acknowledge this concern and will expand the predictor selection section.
  We tested a broader initial set including topographic wetness index (TWI), curvature, drainage density, and NDVI, but these showed Pearson correlations |r| > 0.75 with the retained variables in the correlation and PCA screening, and were excluded to avoid multicollinearity (Reichenbach et al., 2018, Earth-Science Reviews, 180:60–91). We will report the full correlation matrix as supplementary material.
  
  Response to Additional Comments and Minor Comments
  We thank the reviewer for these constructive observations, which we consider will substantially improve the clarity and rigour of the manuscript. We commit to addressing all of them in the revised version.
  
  Citation: https://doi.org/10.5194/egusphere-2026-1749-AC1
- CC1:
  'Reply on RC1. i have a doubt in M3 model . code suggests it has positive Morans value but in paper it is negative.', Basit Ahad Raina, 09 Jun 2026
  
  i have a doubt in M3 model . code suggests it has positive Morans value but in paper it is negative.
  
  Citation: https://doi.org/10.5194/egusphere-2026-1749-CC1
  - AC3: 'Reply on CC1', Edier Vicente Aristizábal Giraldo, 04 Sep 2026
    
    Thanks for your observation. It has been adjusted in the final version.
    
    Citation: https://doi.org/10.5194/egusphere-2026-1749-AC3
RC2:
'Comment on egusphere-2026-1749', Haoyuan Hong, 06 Aug 2026

Review report

Thank you for submitting your manuscript on “The dominant role of latent spatial structure in landslide susceptibility”. The research question is not clear, the idea and method are a little interesting. There are some texts should be improved such as introduction, validation and discussion. Therefore, I give some suggestion and question which I hope useful to the author. My decision is MAJOR REVISION.

1. Introduction

It was difficult to see the justification for the need of this research. The literature review is poor. The paper needs to clearly state what are the problems with the existing works (these types of approaches) and what problem(s) this particularly paper was going to address. Without this clearly problem statement readers would have difficulty to see the merit of this paper. The author only lists some references, I did not find the problem with the exist method. The problem of the existing method is not clear. The author should show us deep analysis about the gap between existing method. Then the research question should be described clearly.

2. Methodology and data

1) The key idea of the method should be added.

2) The flowchart of the method should be added.

3) What is the key novelty of the method should be added.

4) Which sampling method have been applied to generate landslide data and non-landslide data? Did the reliability of the data considered in the model.

5) What is the ratio between train data and test data, why?
6) Which method have been applied to validate the proposed method, is there field survey?

7) What is the detail information about landslide inventory map, such as loss, type, area, some photo of landslide should be added, the origin data format of landside is polygon-based or point based.

3. Results and discussions

There were very few discussions of previous studies. The author should pay more attention or deeper analysis about the effect caused by the parameter of the model.

Other comment:

Figure: The resolution of figure should be improved.

Reference: There are some latest articles should be updated.

Citation: https://doi.org/10.5194/egusphere-2026-1749-RC2
- AC2: 'Reply on RC2', Edier Vicente Aristizábal Giraldo, 04 Sep 2026
  
  We thank the reviewer for the constructive comments, which we address point by point below. Changes made to the manuscript in response to this reviewer are highlighted in green in the revised manuscript.
  
  1. Introduction — research gap and justification
  
  Reviewer comment: "It was difficult to see the justification for the need of this research. The literature review is poor. The paper needs to clearly state what are the problems with the existing works (these types of approaches) and what problem(s) this particularly paper was going to address. Without this clearly problem statement readers would have difficulty to see the merit of this paper. The author only lists some references, I did not find the problem with the exist method. The problem of the existing method is not clear. The author should show us deep analysis about the gap between existing method. Then the research question should be described clearly."
  
  Response: We agree the research gap could be stated more explicitly. We have added a short paragraph at the end of the Introduction that directly signposts: (i) the problem with existing models (the spatial independence assumption inflates predictor significance and produces overconfident estimates, as detailed earlier in the Introduction); (ii) the specific gap addressed here (hierarchical spatial autoregressive models jointly handle heterogeneity and dependence but have never been applied to landslide susceptibility, unlike existing landslide-specific spatial corrections such as eigenvector filtering or autologistic regression, which address dependence alone); and (iii) the explicit research question and objectives, retained in the following paragraph. See the revised Introduction (Sect. 1).
  
  2. Methodology and data
  
  2.1) The key idea of the method should be added.
  
  Response: We have added a brief, plain-language summary of the modeling approach at the start of Section 3 (Data and methods), before the technical/equation-level description, so that readers can grasp the overall strategy without first parsing the GLM/CAR formalism.
  
  2.2) The flowchart of the method should be added.
  
  Response: We have added a workflow diagram as the new Figure 1, immediately after the plain-language summary at the start of Section 3, showing the full pipeline from data inputs (landslide inventory, predictors, catchments) through the model progression (M1-M5), diagnostics, model selection, and the resulting susceptibility map. All subsequent figures have been renumbered accordingly (former Fig. 1-7 are now Fig. 2-8).
  
  2.3) What is the key novelty of the method should be added.
  
  Response: We have added an explicit novelty statement in the Introduction, immediately after the gap statement, distinguishing why this study is needed (the gap) from what is specifically new about it (the novelty): (i) hierarchical spatial autoregressive (CAR-based) models, which jointly represent regional heterogeneity and local dependence in a single Bayesian framework, applied here for the first time to landslide susceptibility; and (ii) a quantitative diagnostic protocol for testing whether predictor-response relationships in a given susceptibility model are confounded by unmodeled spatial structure.
  
  2.4) Which sampling method have been applied to generate landslide data and non-landslide data? Did the reliability of the data considered in the model.
  
  Response: Our approach does not use a landslide/non-landslide (presence/absence) sampling design, since the response variable is a count, the number of landslides per catchment (526 catchments), not a binary classification target. This Poisson count formulation does not require selecting non-landslide sample locations or a sampling ratio: every catchment, whether it contains zero or many landslides, contributes with its full geometry as the exposure unit (via the offset for catchment area). Regarding data reliability, we have added a paragraph in Section 3 (in response to Reviewer 1) explicitly discussing inventory limitations: it was digitized by a single interpreter without independent cross-validation, and high-resolution imagery availability is uneven across the 53-year record and the three basins — both flagged as potential sources of detection bias.
  
  2.5) What is the ratio between train data and test data, why?
  
  Response: The five models were originally compared using in-sample Bayesian goodness-of-fit criteria (DIC, WAIC), which is standard practice for hierarchical spatial model comparison, rather than a manual train/test split. In direct response to this comment, we additionally implemented a spatially-blocked 5-fold cross-validation (catchments grouped into five geographically compact folds via k-means on catchment centroids, following Brenning, 2005, NHESS, rather than randomly, to avoid leakage from spatial autocorrelation between neighboring catchments), with results reported in the new Table S3 (Supplement) and discussed in a new paragraph in Section 5.1. The result is important and is now reported transparently: while M5 (Leroux) remains the best-fitting model in-sample, it and the ICAR model (M3) become numerically unstable when an entire spatial neighborhood is held out, because their predictions rely almost entirely on local spatial smoothing from now-absent neighbors. The non-spatial models (M1, M2) and the BYM model (M4) predict more stably out-of-sample under this harsh test. We discuss this explicitly as a limitation of strongly CAR-dominated models in the revised Discussion.
  
  2.6) Which method have been applied to validate the proposed method, is there field survey?
  
  Response: Our study does not include a field survey component; the landslide inventory was compiled through remote-sensing interpretation, and model validation relies on statistical diagnostics rather than field-based ground-truthing. We validate the models internally via: (1) information-theoretic goodness-of-fit criteria (DIC, WAIC); (2) residual spatial autocorrelation (Moran's I); (3) posterior credible intervals for all fixed effects (Table S2); and (4) spatially-blocked k-fold cross-validation for out-of-sample predictive performance (response to comment 2.5, Table S3). We acknowledge that field verification of a subset of mapped landslides and of the resulting susceptibility zones would further strengthen confidence in the results, and flag this as a priority for future work, consistent with the inventory limitations already discussed in Sect. 3.
  
  2.7) What is the detail information about landslide inventory map, such as loss, type, area, some photo of landslide should be added, the origin data format of landside is polygon-based or point based.
  
  Response: The inventory is point-based (crown-point location), not polygon-based, consistent with common practice for large regional inventories digitized by a single interpreter; we have clarified this in Section 3. This also explains why individual landslide areas were not available for a magnitude-frequency completeness analysis (already noted as a limitation in Sect. 3, in response to Reviewer 1). Landslide type is already described in Section 3 (predominantly shallow translational soil slides and debris flows mobilizing residual soils and saprolite on hillslopes steeper than 25°). We did not collect loss/damage data or field photographs as part of this remotely-sensed inventory; we acknowledge this as a limitation and note it as a priority for a future, dedicated field-based re-mapping effort, consistent with the inventory limitations paragraph already added in Sect. 3.
  
  3. Results and discussions
  
  Reviewer comment: "There were very few discussions of previous studies. The author should pay more attention or deeper analysis about the effect caused by the parameter of the model."
  
  Response: We have addressed the literature-comparison part of this comment already, in response to a closely related comment from Reviewer 1 (Additional Comment 7), by adding a paragraph in Section 5.2 that quantitatively compares our findings (magnitude of residual-autocorrelation reduction, dominance of the structured spatial component) against Lombardo et al. (2019, 2020), Li et al. (2019), and the disease-mapping literature (Lawson, 2018; Morris et al., 2019). For the parameter-effect analysis specifically requested here, we have added a further paragraph in Section 5.2 that interprets the practical (not merely statistical) magnitude of the fixed-effect coefficients: expressed in interpretable raw units, mean slope and elevation have compounding, several-fold effects on expected landslide frequency across the observed range of the study area, whereas the rainfall-day effect is small and, in the more flexible spatial models, not credibly different from zero — indicating that terrain morphometry, not the rainfall regime, dominates the spatial pattern of relative susceptibility in this setting.
  
  Other comments
  
  Figure resolution should be improved.
  
  Response: We thank the reviewer for flagging this. Two panel figures (the landslide inventory map and the covariance-matrix visualization) are at noticeably lower resolution (800x1000 px, 96 DPI) than the rest of the figures (2400-5800 px, 300-500 DPI). The corresponding author will re-export these two panels at higher resolution before final submission.
  
  Reference: There are some latest articles should be updated.
  
  Response: We have added several recent and directly relevant references throughout the manuscript. In the Introduction, we now cite a 2025 study applying the inlabru/INLA framework to landslide susceptibility via spatial point processes (Suen et al., 2025), reinforcing that INLA-based spatial modeling remains an active, current research direction. In the Discussion, we cite a 2025 NHESS brief communication on jointly visualizing susceptibility and its uncertainty (Schlögl et al., 2025), which we relate to our own credible-interval-based uncertainty quantification and identify as a natural extension of this work. Additional recent references were also added throughout in response to Reviewer 1's comments (e.g., Loche et al., 2022; Li et al., 2019).
  
  Citation: https://doi.org/10.5194/egusphere-2026-1749-AC2

Edier Aristizábal, Luigi Lombardo, and Oliver Korup

Data sets

database Edier Aristizabal https://github.com/edieraristizabal/PAPER_BHGLM

Edier Aristizábal, Luigi Lombardo, and Oliver Korup

Viewed

Total article views: 372 (including HTML, PDF, and XML)

HTML	PDF	XML	Total	BibTeX	EndNote
253	94	25	372	23	24

HTML: 253
PDF: 94
XML: 25
Total: 372
BibTeX: 23
EndNote: 24

Views and downloads (calculated since 18 May 2026)

Month	HTML	PDF	XML	Total
May 2026	134	42	8	184
Jun 2026	27	8	5	40
Jul 2026	46	20	8	74
Aug 2026	34	19	3	56
Sep 2026	12	5	1	18

Cumulative views and downloads (calculated since 18 May 2026)

Month	HTML	PDF	XML	Total
May 2026	134	42	8	184
Jun 2026	27	8	5	40
Jul 2026	46	20	8	74
Aug 2026	34	19	3	56
Sep 2026	12	5	1	18

Viewed (geographical distribution)

Total article views: 329 (including HTML, PDF, and XML) Thereof 329 with geography defined and 0 with unknown origin.

Country	#	Views	%

Latest update: 04 Sep 2026

Short summary

Traditional models for landslide susceptibility often ignore spatially-dependent unobserved factors. This study improves accuracy in the Colombian Andes by incorporating spatial dependence and landscape variety, which are natural in environmental data. By accounting for these hidden physical connections, we found that landslide occurrences are more clustered than once thought.


Total:	0
HTML:	0
PDF:	0
XML:	0