# XAI Model for Accurate and Interpretable Landslide Susceptibility

Geophysical Research Letters

Department of Chemistry and Biochemistry, University of California Los
Angeles, 607 Charles E. Young Dr. East, Los Angeles, CA 90095-1569, USA.
Department of Earth, Planetary, and Space Sciences, University of California,
Los Angeles, 595 Charles E. Young Dr. East, Los Angeles, CA 90095-1567
Department of Bioengineering, University of California Los Angeles, 607
Charles E. Young Dr. East, Los Angeles, CA 90095-1569, USA. California
NanoSystems Institute, University of California Los Angeles, 607 Charles E.
Young Dr. East, Los Angeles, CA 90095-1569, USA. Current address: Biomedical
Imaging Research Institution, Cedars-Sinai Medical Center, 116 N Robertson
Blvd, Los Angeles, CA 90048, USA.

∗ co-first authors

S.G. Moon, L.-S<EMAIL_ADDRESS><EMAIL_ADDRESS>
Abstract. Landslides are notoriously difficult to predict. Deep neural
networks (DNNs) models are more accurate than statistical models. However,
they are uninterpretable, making it difficult to extract mechanistic
information about landslide controls in the modeled region. We developed an
explainable AI (XAI) model to assess landslide susceptibility that is
computationally simple and features high accuracy. We validated it on three
different regions of eastern Himalaya that are highly susceptible to
landslides. SNNs are computationally much simpler than DNNs, yet achieve
similar performance while offering insights regarding the relative importance
of landslide control factors in each region. Our analysis highlighted the
importance of: 1) the product of slope and precipitation rate and 2)
topographic aspects that contribute to high susceptibility in landslide areas.
These identified controls suggest that strong slope-climate couplings, along
with microclimates, play more dominant roles in eastern Himalayan landslides.
The model outperforms physically-based stability and statistical models.

##  1 Introduction

Landslides are a major natural hazard that cause billions of dollars in direct
damages and thousands of deaths globally each year [Petley (2012), Froude
Petley (2018)]. Landslides can also cause various secondary hazards, such as
damming and flooding, which often leave a region prone to subsequent damage
following the initial event. [R. Huang Fan (2013)] Additionally, landslide
debris may cause instability by perturbing river sedimentation and disrupting
ecosystems [Fan . (2019), R. Huang Fan (2013)]. As landslide hazards are
expected to increase due to climate change, scientists have sought to more
accurately measure landslide susceptibility [Tien Bui . (2012), Tien Bui .
(2019), Phong . (2019), Dikshit . (2020), Kirschbaum . (2020), Stanley
Kirschbaum (2017)], an estimate of the probability that a landslide may occur
in a specific area, with the goal of mitigating its impact on the economy,
safety, and local ecosystems.

Landslide occurrences are affected by various factors including physical
attributes of the terrain, such as topographic slopes and drainage areas, and
material properties such as the density and strength of soil and bedrock
[Dietrich . (1995), Montgomery Dietrich (1994), Montgomery . (1998), Radbruch-
Hall (1982)]. Also, environmental conditions such as climate, hydrology,
ecology, and ground motion due to earthquakes play an important role in
causing slope instability [Guzzetti . (1999), Baum . (2002), Meunier .
(20081)]. Landslide susceptibility is calculated from these various
controlling factors either through physically based models [Montgomery
Dietrich (1994), Montgomery . (1998), Baum . (2002), Baum . (2010)], data-
driven approaches utilizing statistical analysis [Lee Sambath (2006), Regmi .
(2014)], or machine learning techniques (ML), including support vector
machines and deep neural networks [Tien Bui . (2019), Van Dao . (2020), Tien
Bui . (2020), Conforti . (2014), Gómez Kavzoglu (2005), Lee . (2004)].

While substantial work has been devoted to assessing susceptibility, each
model has shortcomings. Physically- or mechanistically-based approaches, built
on limit equilibrium between driving and resisting forces, have been widely
applied to assess slope stability over a region [Montgomery Dietrich (1994),
Dietrich . (1995), Montgomery . (1998), Dietrich . (2001)]. However,
mechanistic models have limitations, including a limited number of variables,
simplified assumptions of certain environmental conditions (e.g., antecedent
moistures, bedrock structure) and landslide geometry, or the high cost of
geotechnical exploration [Guzzetti . (1999)]. Alternatively, data-driven
approaches, including statistical and ML methods, can handle a large number of
controls to assess susceptibility. Statistical methods such as logistic
regression and likelihood ratios [Lee Sambath (2006), Regmi . (2014),
Reichenbach . (2018)] can utilize a multitude of landslide controls as inputs.
Scientists using these data-driven approaches have obtained a measurable
degree of success in determining areas susceptible to landslides [Lee Sambath
(2006), Regmi . (2014), Tien Bui . (2019)]. However, these data-driven models
rely on the expert’s choice, preconditions, and classifications of input
variables. The outcome of these models’ results, the landslide susceptibility
map, does not decouple individual feature contributions to landslide
susceptibility nor account for their interdependencies due to the limited
computational capabilities in conventional approaches [Reichenbach . (2018)].

Machine learning approaches such as fuzzy logic algorithms, support vector
machines, and DNNs have been applied to landslide studies for mapping
landslide susceptibility [Pradhan (2013), Tien Bui . (2020), Gómez Kavzoglu
(2005)]. DNNs have achieved superior performance compared to both statistical
methods and other ML approaches due to their use of nonlinearities, complex
interdependencies of interlayer connections, as well as internal
representations of data [Conforti . (2014), Gómez Kavzoglu (2005), Tien Bui .
(2020), Van Dao . (2020), Rudin (2019), Gunning . (2019), Adadi Berrada
(2018)]. However, the “black box” nature of DNNs has been a major hurdle for
their adoption in practice and research, making it difficult for experts to
understand and trust the outcome. With DNNs it is nearly impossible to
determine the exact relation between individual inputs and outputs [Rudin
(2019), Gunning . (2019), Adadi Berrada (2018)]. Lack of interpretability is a
weakness of DNNs and a fundamental drawback for high-stakes applications such
as landslide mitigation where decisions impact lives and result in untold
costs of insurance and reconstruction [Cui . (2019), Froude Petley (2018), R.
Huang Fan (2013)].

The DNN’s lack of interpretability has prompted Defense Advanced Research
Projects Agency (DARPA)’s third wave of AI call in 2017 and the European
Union’s 2018 General Data Protection Regulation, which grants a “right to an
explanation,” for algorithmic decisions that are made [European Commission
(2020)]. Next-generation AI systems refers to the so-called explainable or
interpretable AI (XAI) models. The latter must be able to construct
explanatory models for classes of real-world phenomena that can be
communicated to humans [Adadi Berrada (2018)]. Various XAI categories have
since been defined in the literature based on factors such as application and
methodology, where each category is further divided into subclasses [Li .
(2020)]. Although the use of XAI in research is expanding, existing approaches
aimed at explaining black box models exhibit a trade-off between accuracy and
interpretability, resulting in a large gap in performance. Recently, Rudin
[Rudin (2019)] showed that with proper feature engineering, and a shift from
explaining existing black box models to creating methods with inherently
interpretable models, the trade-off between accuracy and interpretability can
be circumvented.

There is an urgent need for XAI methods in the analysis of geohazards [Bergen
. (2019), Toms . (2020), Dikshit . (2020)]. At this time, none exist
specifically for the modeling of landslide hazards. Thus, we developed a novel
end-to-end method that bridges the gap between explainability and accuracy.
Our XAI method, termed superposable neural network (SNN), is data-driven, in
the sense that explanations are generated purely from the data without
external information, hardwired physical relationship, or preconditioning. It
is a hybrid between model extraction methods and feature-based methods. Model
extraction methods aim to train an explainable “student” model to mimic the
behavior of a “teacher” model, which is taken here to be a high performing
pre-trained black box model [Li . (2020)]. On the other hand, feature-based
methods aim to analyze the influence or the importance of input features [Li .
(2020)].

##  2 Data and Methods

###  2.1 Superposable neural networks

Our XAI model is an additive neural network architecture enforcing no
interconnections between inputs. The lack of interconnections between features
is what enables explainability. Additive models have been gaining increasing
interest in XAI [Agarwal . (2020), Hastie Tibshirani (1990), Friedman (2001)],
where separate neural networks each specialized in a single feature are
combined to optimize a common outcome. Such models enable determining the
contribution of each input feature to the final outcome. However, state-of-
the-art algorithms rely on gradient descent-based optimizers to obtain model
parameters, making them notoriously difficult to train jointly and expand
beyond a limited number of features. Herein we introduce a new training
approach that enables the separate training of individual neural networks. Our
approach utilizes several state-of-the-art ML techniques (multi stage
training, knowledge distillation, second order optimization [Youssef . (2015),
Bouchard Youssef (2018), Youssef . (20181), Yu Wilamowski (2011), Hinton .
(2015)]) to deliver a model that is optimal in terms of performance and
remarkably simple in terms of architecture.

In the superposable neural network (SNN; Figs. 1 and S1), instead of
interdependencies being embedded in layers of network connections,
interdependencies are explicitly created as a product function of multiple
original input features, which we refer to as composite features. An automated
feature selection method finds interdependencies between features prior to
training the network by isolating composite features that contribute to the
desired outcome. Contributing composite features are explicitly added as
independent inputs to the model, while non-contributing composite features are
discarded (see Supporting Information).

Composite features are categorized by the number of product operations
involved. For example, given a problem with $n$ original input features
$x_{1},x_{2},...x_{n}$, we can generate a set of $M\geq n$ composite features
$\chi_{1},\chi_{2},...\chi_{M}$, where Level-1 features are merely the single
original features (first-degree monomials such as $x_{i}$). Level-2 features
are composite features defined as the product of two first-degree monomials.
As an example, we may form the product $x_{1}*x_{2}$ (second-degree monomial),
where the monomials $x_{1}$ and $x_{2}$ are Level-1 features. Level-3 features
are composite features consisting of a product of three first-degree
monomials, such as $x_{1}*x_{2}*x_{3}$ or $x_{1}*x_{2}^{2}$, and so on,
resulting in third-degree monomials. Composite features are restricted to
functions that cannot be derived from another function by elementary algebraic
transformations. For example, $x_{1}^{2}*x_{2}^{2}$ and $2*x_{1}*x_{2}$ are
not permitted since they can be derived from $x_{1}*x_{2}$ by elementary
operations (namely, by squaring and scaling, respectively). Using composite
features, SNNs can approximate any continuous function for inputs within a
specific range as a polynomial expansion to any desired precision, giving SNNs
the ability to retain a level of accuracy on par with state-of-the-art DNNs.
While this architecture results in outputs where the contributions from each
input feature can be separated, we are left with the non-trivial question of
how can one train such a deceptively simple SNN to yield high performance on
par with state-of-the-art deep nets. The latter are orders of magnitude more
complex in terms of number of parameters and connectivity, making the training
of SNNs appear to be impossible.

The SNN is represented mathematically by the function (Eq. 2):

$\displaystyle
S_{t}(\\{\chi_{j}\\})=\sum_{j}\Bigl{(}\sum_{k}w_{j,k}e^{-(a_{j,k}\chi_{j}+b_{j,k})^{2}}+c_{j}\Bigr{)}.$
(1)

It contains only two hidden layers of neurons with radial basis activation
functions in the first layer and linear activation functions in the second
layer. The choice of radial basis activation functions allows the user to
minimize the number of neurons in the model, which is critically important for
maximizing the efficiency of our method. Each input $\chi_{j}$ is exclusively
connected to a group of neurons to form an independent function
$S_{j}=\sum_{k}w_{j,k}e^{-(a_{j,k}\chi_{j}+b_{j,k})^{2}}+c_{j}$ and the SNN
output $\mbox{$S_{t}$}=\sum_{j}S_{j}$ is the sum of all independent functions,
where $j$ = 1 : number of features ($m$), $k$ = 1 : $v$, the number of neurons
per feature, and $\chi_{j}$ is the $j\textsuperscript{th}$ composite feature.
In addition to determining the features and interdependencies between features
that contribute to the outcome, the SNN architecture enables the
quantification of the exact contribution of each feature and feature
interdependency to the outcome.

Using second-order optimization in lieu of the commonly used gradient based
first-order optimization is a key factor in minimizing model complexity and
maximizing accuracy in our method. Second-order optimization is more robust
with respect to hyperparameter tuning, and offers better accuracy, convergence
speed, and efficiency [Ozyildirim Kiran (2020), Le . (2011), Battiti (1992),
H. Tan Lim (2019), Montavon . (2012), Yu Wilamowski (2011), Wilamowski Yu
(2010)]. However, second-order optimization is not commonly used due to the
high computational and memory storage costs involved. In previous studies, we
have tackled this problem by developing the multistage trained (MST) DNN
[Youssef . (2015), Bouchard Youssef (2018), Youssef . (20181)]. MST is a
special type of deep net consisting of smaller subunits that are trained in
hierarchical fashion using second order optimizers. MSTs employ a hierarchical
structure of multilayer perceptrons, each trained in successive stages,
enabling a better solution to be obtained, requiring only modest compute
resources [Youssef . (2015), Bouchard Youssef (2018), Youssef . (20181)].

In this work, a second-order MST DNN [Youssef . (2015), Bouchard Youssef
(2018), Youssef . (20181)] is used as the teacher model of the SNN. Two stages
of knowledge distillation implemented in a step-wise manner are used to
further facilitate the optimization of the highly constrained SNN architecture
in a way that maximizes accuracy while minimizing the number of neurons for
optimal explainability. The robustness characteristic of second-order
optimization is key to enabling the automation of our feature selection
method, which involves training thousands of DNN models. This would not have
been possible with first-order optimization due to the high sensitivity to
hyperparameters and the manual effort involved in tuning. SNNs based on
knowledge distillation and second-order optimization techniques aim to deliver
an explainable neural network model without compromising accuracy (see details
in Supporting Information).

###  2.2 Study Area

Asia holds the majority of human losses due to landslides globally, with a
high concentration in the Himalayan Arc [Petley (2012), Froude Petley (2018)].
Previous studies in the Himalaya regions presented various controlling factors
for landslides, which include steep slopes, extreme precipitation events,
flooding, and frequent earthquakes [Larsen Montgomery (2012), Coudurier-
Curveur . (2020), Bookhagen Burbank (2010), Chowdhuri . (2021), Kent Dasgupta
(2004)]. The eastern end of Himalaya has numerous landslides, which may be
affected by extreme precipitation events or earthquakes from active tectonics
(Figs. 2, S2). This area exhibits a dramatic precipitation gradient due to
moisture originating from the Bay of Bengal in the south [Bookhagen Burbank
(2010), Barros . (2004), Yang . (2018)] (Fig. S3). Mean annual precipitation
rates (_MAP_) vary from $\scriptstyle\sim$7000 mm/yr in the rangefront to
$\scriptstyle\sim$200 mm/yr in the hinterland [Bookhagen Burbank (2010)]. The
number of extreme rainfall events (_NEE_) calculated as the number of days
that exceed the 90th percentile of daily rainfall rates, reaches
$\scriptstyle\sim$13 and $\scriptstyle\sim$2 events/yr in the rangefront and
hinterland, respectively [Bookhagen Burbank (2010)]. Both daily and mean
annual precipitation rates are calculated based on 3-h measurements from the
Tropical Rainfall Measuring Mission (TRMM) 2B31 over 12 years (January 1998 to
December 2009), with a spatial resolution of $\scriptstyle\sim$2.6 km
[Bookhagen Burbank (2010)]. Orographic patterns of precipitation rates and
variability are also observed in the 57-yr Asian Precipitation–Highly Resolved
Observational Data Integration Towards Evaluation of Water Resources project
(APHRODITE) [Yatagai . (2012)]. In addition, this area has consistently steep
slopes (>0.59), from the rangefront, where Holocene Himalayan shortening is
concentrated near and along the Main Frontal Thrust, into the hinterland,
which is affected by deglaciations from the last glacial maximum [Burgess .
(2012), Haproff . (2019), Haproff . (2020), Salvi . (2020)]. Furthermore, the
area is prone to active seismicity. The 1950 Ms 8.6 Assam earthquake, one of
the largest earthquakes in the Himalayan range, struck near the Namche Barwa
region [Ben-Menahem . (1974)]. Additionally, the eastern end of Himalaya has
experienced $>$450 earthquakes with M$>$4 since 1973 according to Incorporated
Research Institutions for Seismology data archive (www.iris.edu, accessed on
10/01/2020).

We selected three regions (N-S, NW-SE, and E-W regions) within the eastern end
of Himalaya that have varying ranges of landslide controls to test the
performance of the SNN model (Figs. S2, S4). Both the N-S and E-W regions
extend from the active rangefront to the hinterland, from north to south and
east to west, respectively. The NW-SE region mostly lies along the active
rangefront. Testing the SNN over three regions will allow us to examine the
following: 1) whether there are universal or distinctly different controls of
landslides; and 2) whether SNN-determined independent functions of feature
contributions to susceptibility, $\\{S_{j}\\}$, are similar across these three
regions.

###  2.3 Landslide Inventory

We generated a landslide inventory of the eastern end of Himalaya by combining
manual delineation by detecting vegetation changes with a semi-automatic
detection algorithm [Ghorbanzadeh . (2019), Prakash . (2020)]. We used 30 m
resolution Landsat 8 imagery from November 2017 with bands 2, 3, 4, 5, and 7.
[United States Geological Survey EarthExplorer (2020)] These satellite images
were used to generate natural and false color imagery to best identify
landslides. High degrees of vegetation in the area allow for clear delineation
of landslides suitable for mapping. Then, landslide polygon detection and
delineation were mainly based on Landsat spectral signatures, with
consideration for topographic and climatic information such as slope, relief,
and _MAP_ in an automatic detection algorithm [Krizhevsky . (2012), Lawrence .
(1997), Ghorbanzadeh . (2019), Prakash . (2020)]. Next, the output from the
automatic detection algorithm was visually compared with Landsat 8 imagery,
high-resolution Google Earth images, and a 4-band Planetscope Scene with a 3 m
resolution. We manually corrected landslides from the automatic detection
method to resolve common problems such as the inaccurate representation of
landslide areas because of aggregated features and small detected features
that are not related to landslides [Marc Hovius (2015), Parker . (2011)]. We
used the manually corrected, automatically mapped landslides as our final
landslide inventory (details in Supporting Information).

Our landslide inventory covers our overall study area of 4.19$\times$109 m2
and has 2,289 landslides whose areas range from 900 to 1.96$\times$106 m2.
Landslide densities calculated over a 2.25 km2 window are high in the
rangefront (maximum of 0.121) and low in the hinterland (maximum of 0.039)
(Fig. 2). The area-frequency distribution of landslides is similar to a
previous study from a nearby region [Larsen Montgomery (2012)] (Fig. S3). They
showed that a transition from soil to bedrock landslides occur at a landslide
area of $\scriptstyle\sim$100,000 m2 [Larsen . (2010), Larsen Montgomery
(2012)]. In our landslide inventory, $<$1% of landslides in numbers and $<$20%
of total landslide area are greater than 100,000 m2 (Table S1). Thus, most of
the mapped landslides are likely soil landslides. The spatial distributions
and extents of landslides from our inventory are shown in Fig. 1a, Figs. S5a
and S6a, for the N-S, NW-SE, and E-W regions, respectively. Landslide density
is calculated within a 2.25 km2 window, which is greater than the largest
landslide size (1.95 km2) (Fig. 2).

###  2.4 Feature Descriptions

We quantified the spatial distribution of 15 topographic, climatic, and
geologic controls and used them as input features for SNN. The 15 single
features include topographic controls of aspect (_Asp_), mean curvature
(_Curv_ M), planform curvature, profile curvature, total curvature, distance
to channel (_Dist_ C), drainage area, elevation (_Elev_), local relief
calculated as an elevation range within a 2.5 km radius circular window
(_Relief_), and slope, climatic or hydrologic controls of discharge, mean
annual precipitation (_MAP_), and number of extreme events (_NEE_), and
geologic controls of distance to lithologic boundaries (i.e., mostly
faults)(_Dist_ F) and distance to the Main Frontal Thrust and suture zone
(_Dist_ ${}_{\emph{MFT}}$). These features were selected from literatures that
examine landslide occurrences in the Himalayan region [Devkota . (2013), Regmi
. (2014), Mandal Mandal (2018), Chowdhuri . (2021)]. We mostly used features
directly measured through satellite data such as the 90 m digital elevation
model from the Shuttle Radar Topography Mission (SRTM) [United States
Geological Survey EarthExplorer (2020)], rainfall rates and variability from
TRMM [Bookhagen Burbank (2010)], and published geologic maps [Haproff .
(2019), Taylor Yin (2009)]. Utilizing the open-source satellite data with a
long-term historic archive allows our approach to be easily implemented in
other regions (e.g., Himalayan Arc) with limited accessibility, high landslide
potential, and long landslide histories [Zhu . (2019), Petley (2012), Froude
Petley (2018), Kirschbaum . (2020)]. The details of our data sources and
methods of calculation as well as the uncertainties and limitations of our
study are explained in the Supporting Information (Fig. S3, Table S2).

###  2.5 Assessments of Model Performances

A threshold value of $S_{t}$ can be used as a binary classifier to predict
landslides and compare them with observed landslides from our inventory. We
selected a threshold susceptibility corresponding to the closest point to a
perfect classifying model of a 100% true positive rate and 0% false positive
rate in a receiver operating characteristic (ROC) curve. Areas with $S_{t}$
greater and lower than this threshold are classified as landslide (_ld_) and
non-landslide (_nld_) areas, respectively, in model. We used the threshold
$S_{t}$ values of 0.406, 0.519, and 0.572, for the N-S, NW-SE, and E-W
regions, respectively.

We assessed SNN performance based on various metrics including area under the
receiver operating characteristic curve (AUROC) and success rate curve. In
addition, we calculated the statistical measures of accuracy, sensitivity
(probability of detection, POD), specificity (probability of false detection,
POFD), and POD-POFD. The performance and explainability of the SNN were
compared with results from 1) a physically-based model (SHALSTAB) suitable for
soil landslides, which are the dominant landslide type in the region
[Montgomery Dietrich (1994), Dietrich . (2001), Moon . (2011), Larsen
Montgomery (2012)]; 2) the most commonly used statistical models, logistic
regression and likelihood ratios [Lee (2005), Akgun (2012), Reichenbach .
(2018)]; and 3) a conventional deep net trained using the adaptive moment
estimation (Adam) optimizer [Kingma Ba (2014)]. We also calculated the 95%
confidence interval of mean AUROC from the statistical and neural network
model outputs based on a 10-fold cross validation (details in Supplemental
Information).

###  2.6 Identifying Landslide Controls and Quantifying Their Contributions

We used our SNN outputs to better understand the relative importance of
controlling factors and their functional relationships with landslide
susceptibility. By utilizing the exact contribution of $S_{j}$ from each
feature at the pixel scale, we quantified the relative importance of landslide
controls vary in different localities. We identified primary controls of
landslides that induce large differences in average $S_{j}$ between (_ld_) and
_nld_ areas (i.e., $\Delta\bar{S}_{j}$ = $\bar{S}\textsubscript{\emph{j,ls}}$
\- $\bar{S}\textsubscript{\emph{j,nls}}$)) within a given area (e.g., 2.25
km2) (Fig. 2, Figs. S5, S6). Additionally, we calculated $\Delta\bar{S}_{j}$
using observed landslide and non-landslide areas from our inventory to compare
with those from our SNN-predicted landslides. Both results show similar
rankings and magnitudes of $\Delta\bar{S}_{j}$ in our studied regions (Fig.
2d-e, Figs. S5d, S6d). For comparison, we also identified primary controls of
landslides and their relative contributions from SNN-Level1 and the weights
determined by the logistic regression model (Figs. S7, S8). Lastly, we
examined the independent functions $\\{S_{j}\\}$ to better understand the
controls of individual features on landslide susceptibility. The details of
our methods are explained in the Supporting Information.

##  3 Results

The SNN model performance is comparable to that of second-order-optimized
DNNs, while providing a significant improvement over commonly used physically-
based model, statistical models, and first-order-optimized DNNs (Tables S3).
AUROCs of Level-1 and Level-2 SNNs are 0.864 and 0.896 (N-S), 0.826 and 0.856
(NW-SE), and 0.878 and 0.919 (E-W), with Level-2 SNNs capturing over 98% of
the teacher model (MST) performance across all three study regions.
Additionally, we tested Level-3 SNN in one of the regions (NW-SE) and found it
to perform marginally better ($\scriptstyle\sim$1%) than Level-2 SNN. Given
that the performance of Level-2 already exceeds 98% of the teacher model’s
performance, we adopt the Level-2 SNN for this study. The Level-2 SNN is
optimal in the sense that it provides high accuracy (comparable to deep nets)
and relatively simple model complexity (hereafter, SNN refers to Level-2 SNN).
Our SNN shows a statistically significant improvement over a conventional DNN
(AUROC = 0.837 (N-S), 0.762 (NW-SE), 0.835 (E-W)). The SNN achieved 19-23%
improvement in AUROC over the top performing single original features of _MAP_
or slope (AUROC = 0.756 (N-S) or 0.696 (NW-SE) and 0.760 (E-W), respectively),
19-26% improvement over a physically-based model (SHALSTAB) (AUROC = 0.756
(N-S), 0.694 (NW-SE), 0.732 (E-W)), 3-10% improvement over logistic regression
(AUROC = 0.869 (N-S), 0.795 (NW-SE), 0.880 (E-W)) and likelihood ratios (AUROC
= 0.841 (N-S), 0.792 (NW-SE), 0.835 (E-W)) in our three study regions. The 95%
confidence intervals of the mean AUROC of the SNN lie above and do not overlap
with those of the statistical models and DNN (Tables S4). In addition, other
performance metrics such as accuracy, POD, POFD, and POD-POFD from the SNN are
improved over other methods as well (Tables S5).

We identify composite features as locally important landslide controls on our
study area. Namely, the product of slope times _NEE_ or _MAP_ , and _Asp_
(Fig. 2c, d). The primary control of _MAP*Slope_ and _NEE*Slope_ is consistent
among our three study regions in the eastern end of Himalaya (Fig. 2c, Figs.
S5c, S6c) despite differences in the spatial distribution and magnitude of
precipitation and proximity to a major fault with a history of earthquakes
(Fig. S4). The overall importance of slope and climatic features identified by
the SNN are similar to those inferred from the weights determined by logistic
regression. However, the relative and absolute importance of features
decoupled from other features by the SNN are different from those determined
by logistic regression (Fig. S7). We note that the model performances obtained
here cannot be directly compared to DNN performances cited in the literature
because these models are trained on different regions; this is why independent
validation against competing methods are provided here over the same region.

To assess the relative importance between climate and slope features, we used
the SNN Level-1 results. We first add all single original features’
$\Delta\bar{S}_{j}$ related to climate ($j$ = _Asp, NEE, MAP_) (hereafter,
$\Delta\bar{S}_{\emph{climate}}$) and compared it to that of slope,
$\Delta\bar{S}\textsubscript{\emph{Slope}}$, for a given area. Average ($\pm$
1$\sigma$) $\Delta\bar{S}_{\emph{climate}}$ and
$\Delta\bar{S}\textsubscript{Slope}$ values are 0.044 $\pm$ 0.028 and 0.052
$\pm$ 0.050 (N-S), 0.042 $\pm$ 0.024 and 0.052 $\pm$ 0.040 (NW-SE), and 0.031
$\pm$ 0.023 and 0.059 $\pm$ 0.040 (E-W), respectively. Interestingly,
$|\Delta\bar{S}_{\emph{climate}}|$ exceeds $|\Delta\bar{S}_{\emph{Slope}}|$ in
$\scriptstyle\sim$59%, $\scriptstyle\sim$53%, and $\scriptstyle\sim$38% of
localities, which encompass $\scriptstyle\sim$84%, $\scriptstyle\sim$60%, and
$\scriptstyle\sim$59% of landslide areas in the N-S, NW-SE, and E-W regions,
respectively (Fig. S8). We find that the importance of climate feature
contributions to landslide occurrence can exceed that of slope by
$\scriptstyle\sim$30 times in certain locations and in most localities with
high landslide density.

The SNN-determined independent functions $\\{S_{j}\\}$ show varying
relationships between features, feature interdependencies, and both their
absolute susceptibility contribution (Fig. 3, Figs. S9, S10).
$S\textsubscript{\emph{MAP*Slope}}$, $S\textsubscript{\emph{NEE*Slope}}$, and
$S\textsubscript{\emph{Elev*Slope}}$ generally exhibit steep increases with
feature value, followed by asymptotic behavior for the former two (Fig. 3a,
Figs. S9a, S10a). The nonlinear relationships between landslide susceptibility
and the product of slope and climatic features of _MAP_ and _NEE_ are similar
to the predicted patterns of slope stability controlled by subsurface
saturation (Supporting Information). In addition,
$S\textsubscript{\emph{Asp}}$ shows a peak around $145^{\circ}$, to
$180^{\circ}$, indicating a preference for south-facing slopes, likely due to
moisture from the Bay of Bengal [Bookhagen Burbank (2010)] (Fig. 3b,
Supporting Information). These functional relationships are similar to those
deduced by the statistical method of likelihood ratios ($LR$) that represent
the likelihood of landslide occurrence. However, unlike likelihood ratios,
which assume the same, average likelihood ($LR=1$) for each feature, $S_{j}$
corresponding to $LR=1$ varies depending on a feature’s absolute, decoupled
contribution to landslide susceptibility.

##  4 Discussion

Whereas explainable AI efforts involve a trade-off between accuracy and
interpretability, the new SNN does not compromise accuracy. Given the SNN’s
inherent and unique ability to decouple individual feature contributions and
select feature interdependencies, we can easily isolate the local
contributions from primary controls discovered by the SNN (Fig. 2). Our local
analyses for assessing landslide controls indicate that the contribution of
climate features, such as _NEE_ , _MAP_ , and _Asp_ , to landslide
susceptibility tends to surpass that of slope for a majority of landslide
occurrences in this area. Our findings are different from previous studies in
other Himalayan regions that emphasize the dominance of tectonic and
topographic contributions to erosional processes [Coudurier-Curveur . (2020),
Scherler . (2014)]. Our study highlights a prevalent climatic control on
landslide occurrences in the eastern Himalayan region. Due to the eastward
increasing trends of precipitation rate and variability along the Himalayas,
the eastern end of Himalaya contains one of the largest strike-perpendicular
climatic variations across the steep mountain range [Bookhagen Burbank
(2010)]. This considerable climate gradient in the eastern end of Himalaya
likely impacts landslide susceptibility.

The transparency of our XAI model offers insight into controlling factors of
landslides while offering an accurate quantification of susceptibility.
Specifically, SNNs highlighted important, yet under-appreciated controls of
_NEE*Slope_ , _MAP*Slope_ , and _Asp_ in our study site. The multiplication of
two Level-1 features produces a Level-2 feature that measures the importance
of both features. These controls of the product between slope and
precipitation rates and intensity as well as aspect imply that landslides are
affected by strong slope-climate couplings and aspect-related microclimates.
The nonlinear asymptotic function of $S\textsubscript{\emph{MAP*Slope}}$ and
$S\textsubscript{\emph{NEE*Slope}}$ can be explained by a physical mechanism
of rainfall-induced landslide that induces slope failure due to an increase in
pore-water pressure and subsurface saturation [Iverson (2000)] (details in
Supporting Information). However, this deduced relation is different from
observations of earthquake-induced landslide occurrences that increase
nonlinearly with increasing slope [Meunier . (20082), AL. Huang Montgomery
(2014)]. The identified controls of _MAP_ , _NEE_ , and _Asp_ imply that local
precipitation infiltration affected by the spatial patterns of precipitation
rates as well as moisture availability on steep slopes may be the dominant
contributors to subsurface saturation in the eastern end of Himalaya
(Supporting Information).

Identifying the exact trigger for a landslide requires dense measurements and
historic records of the soil, hydrologic, and climatic conditions (e.g., soil
moisture, antecedent rainfall, rainfall intensity) [Orland . (2020),
Kirschbaum . (2020)], which are often difficult to obtain, especially in rural
mountainous areas with limited accessibility. We have shown that our SNN model
can identify key controls and quantify their potential contributions to
susceptibility, highlighting the essence of strong slope-climate coupled
controls on landslide occurrences. The composite features identified by the
SNN such as _NEE*Slope_ or _MAP*Slope_ are consistent with previous
understanding of landslide mechanisms (Supporting Information). However, they
were not explicitly implemented in previous data-driven statistical models. In
DNNs, such couplings were likely identified, but if that were the case, the
information would have been implicitly contained in the network weights and
not readily available to the user. By incorporating the climatic controls from
_MAP*Slope_ and _Asp_ in the N-S region, the performance of the SNN improved,
increasing AUROC by 3–26% compared to those of statistical or physically-based
models [Montgomery Dietrich (1994), Dietrich . (2001), Lee (2005), Akgun
(2012)] (Supporting Information). This performance enhancement is substantial
and could impact decisions regarding landslide hazards.

Our work presents a significant advance in XAI applications to natural hazards
and circumventing the “black box” nature of common AI models. SNNs provide
quantitative analysis of the controlling factors and further highlight the
important, mechanistic interpretation of landslides. Our new approach based on
tournament ranking and knowledge distillation provides a comprehensive
framework that allows the user to examine numerous composite features,
identify key controls, and more accurately model landslide susceptibility. As
natural perturbations increase due to urban development and climate change,
the SNN may provide a promising, data-driven predictive tool that will enable
communities to tailor plans for hazard mitigation.

We note that our study has limitations. For one, our landslide and original
feature data from the eastern end of Himalaya have relatively coarse
resolutions and are temporally limited (e.g., 30 m resolution Landsat
satellite images from 2017 [United States Geological Survey EarthExplorer
(2020)], 90 m resolution SRTM DEM [United States Geological Survey
EarthExplorer (2020)], and $\scriptstyle\sim$2.6 km resolution TRMM data over
12 years [Bookhagen Burbank (2010)]) (see details in Supporting Information).
Due to the extremely rugged mountains in the Himalayas, the highest available
DEM resolution without extensive data gaps, suitable for regional-scale
landslide susceptability analysis, is 90 m [Kirschbaum . (2020), Stanley
Kirschbaum (2017)]. Also, there are no readily available time-series
precipitation data with a resolution $<$2.6 km in this area. Our study site
doesn’t have detailed information on high-resolution topography, surface
materials (e.g., soil depth, bedrock structures), and climatic and
ecohydrologic conditions (e.g., landslide-triggering storm intensity,
vegetation types). Thus, our results from both physically-based or data-driven
models may be biased due to the inherited uncertainties and limitations of our
input data and may change when sufficient data become available [Leonarduzzi .
(2021)].

Despite data limitations and uncertainties, our SNN analysis of the eastern
end of Himalaya presents an important contribution to landslide hazard
studies. High mountains in Asia hold the majority of human losses due to
landslides globally, with a high concentration in the Himalayan Arc, according
to a global analysis conducted using 2004 - 2016 data [Petley (2012), Froude
Petley (2018)]. Due to the associated high risks, there has been efforts to
model landslide susceptibility in the Himalayan regions based on currently
available data with limited resolutions [Devkota . (2013), Regmi . (2014),
Mandal Mandal (2018), Chowdhuri . (2021)]. Our work aims to capture the
regional-scale, spatial distributions of landslide occurrences, differentiate
controls of landslide hazards, and provide interpretable, essential empirical
relationship between landslide controls and regional-scale landslide
susceptibility. The decoupled, SNN identified functions combined with future
changes in environmental conditions (e.g., extreme precipitation) may provide
a promising tool for assessing the potential landslide hazards in this area
(e.g., [Kirschbaum . (2020)]). Our work is the first application XAI to
landslide hazard, The method is is easily applicable to other locations,
different datasets, and other physical hazards, such as earthquakes and
wildfires. The SNN is remarkably simple consisting of only two hidden layers
yet its performance rivals that of DNNs. Our SNN can be also easily updated
and improved when global, open-source, high-resolution datasets and high-
performance computational resources become more available in the future.

Figure 1: Conventional DNN architecture vs SNN architecture. In a conventional
DNN, features are interconnected and interdependencies are embedded in the
network, making them virtually impossible to separate. In a SNN, features and
feature interdependencies that contribute to the output are found in advance
and explicitly added as independent inputs. The symbols in this figure are
defined in the main text, Equation 2. Details of the SNN architecture are
explained in section 2.1.

Figure 2: Spatial distribution of landslides, susceptibility, and important
controls in the eastern end of Himalaya. Maps are shown for (a) 959 numbers of
mapped landslides in red polygons, (b) total susceptibility at the pixel scale
($S_{t}$) from the Level-2 superposable neural network, and (c) primary
features identified as locally important controls of landslides for the N-S
region of the eastern end of Himalaya. A location of the region is shown in
Fig. S2. The locally important control is identified as the feature with the
largest difference in average contribution ($\Delta\bar{S}_{j}$) between areas
of landslides ($ld$) and non-landslides ($nld$) within a 2.25 km2 window.
Windows of no data in (c) are due to the lack of landslide mapping, shown in
gray dashed outlines, or modeled landslides. Bar charts represent
$\Delta\bar{S}_{j}$ in descending order, and pie charts represent average
$S_{j}$ ($\bar{S}_{j}$) contributions to $ld$ and $nld$ areas for (d) an
example of the local area shown with a white square in (c) and (e) the entire
region. Bars with solid black and gray dashed borders show $\Delta\bar{S}_{j}$
based on modeled and manually mapped landslide inventories, respectively.
Features related to topography, aspect, climate, and geology are shown in
green, pink, blue, and brown or combinations thereof, respectively. Features
outside of those are shown in gray. Mean annual precipitation (_MAP_), aspect
(_Asp_), elevation (_Elev)_ , distances to channel (_Dist C_), all faults
(_Dist F_), and main frontal fault and suture zone (_Dist MFT_), and mean
curvature (_Curv M_). The asterisk * indicates algebraic multiplication of two
features. Information regarding features and variables is provided in
Supporting Information.

Figure 3: Feature contribution to susceptibility ($S_{j}$). Independent
functions of $S_{j}$, identified as primary landslide controls for the N-S
region, are shown for a) mean annual precipitation (_MAP_) * _Slope_ , b)
aspect (_Asp_), and c) elevation (_Elev_) * _Slope_. Likelihood ratios ($LR$)
represent the likelihood of landslide occurrence for a specific range of
feature values, which is calculated as the ratio of the probability of
landslide occurrences to the probability of total area for the given feature
range. $LR$ = 1 and $LR$ $>$ 1 represent the average and above-average
likelihood of landslide occurrence, respectively. $LR$ are shown as short
dashed lines with corresponding right-side $y$-axes. Note that $S_{j}$
corresponding to $LR=1$, shown as a long-dashed gray line, differ between
features because the SNN quantifies the absolute contributions of $S_{j}$
decoupled from other features. Colors are the same as those in Fig. 2c.

## Open Research

For mapping landslides and quantifying landslide controls, we use publically
available satellite data such as 30 m resolution Landsat 8 imagery, the 90 m
digital elevation model from the Shuttle Radar Topography Mission (SRTM)
[United States Geological Survey EarthExplorer (2020)], rainfall rates and
variability from TRMM [Bookhagen Burbank (2010)], and published geologic maps
[Haproff . (2019), Taylor Yin (2009)]. The details of the models are provided
in the Supporting Information.

###### Acknowledgements.

L.-S.B. would like to thank NVIDIA for donating the Quadro P6000 GPU for this
project. S.M. acknowledges the support of NSF EAR-1945431. We thank George
Hilley, Dino Bellugi, Eitan Shelef, and Justin Higa for discussions and Ben
Mirus for constructive review of an earlier version of the manuscript.
Supporting Information for “XAI Model for Accurate and Interpretable Landslide
Susceptibility” Contents

  1. 1.

Text S1. SNN Training Method

  2. 2.

Text S2. Description of Landslide Inventory

  3. 3.

Text S3. Feature Descriptions and Limitations

  4. 4.

Text S4. Explanations for the Controls on Landslide Occurrences

  5. 5.

Text S5. Assessments of Model Performances

  6. 6.

Text S6. Toy model

  7. 7.

Table S1. Description of Landslide Inventory

  8. 8.

Table S2. Description and Ranges of 15 Single Features

  9. 9.

Table S3. Performance metrics for Models and Single Features

  10. 10.

Table S4. Artificial Neural Network and Statistical Model Confidence Intervals

  11. 11.

Table S5. Performance Metrics

  12. 12.

Table S6. Correlation Metrics Between Features (R-value)

  13. 13.

Table S7. Logistic Regression Control Coefficients

  14. 14.

Table S8. Truth table.

  15. 15.

Table S9. Composite features.

  16. 16.

Figure S1. Superposable neural network training flow diagram.

  17. 17.

Figure S2. Elevation map of the eastern end of Himalaya.

  18. 18.

Figure S3. Spatial distribution of 15 features used in the SNN model.

  19. 19.

Figure S4. Comparison between the feature ranges in three regions.

  20. 20.

Figure S5. Spatial distribution of landslides, susceptibility, and important
controls for the NW-SE region.

  21. 21.

Figure S6. Spatial distribution of landslides, susceptibility, and important
controls for the E-W region.

  22. 22.

Figure S7. Bar charts representing $\Delta\bar{Sn}_{j}$ from different
methods.

  23. 23.

Figure S8. Local analysis of climate vs slope susceptibility contributions to
landslide occurrence .

  24. 24.

Figure S9. Primary landslide controls for the NW-SE region are shown.

  25. 25.

Figure S10. Primary landslide controls for the E-W region are shown.

  26. 26.

Figure S11. The relationship among aspect, normalized difference vegetation
index, and ${S}_{Asp}$

  27. 27.

Figure S12. Feature ranking for Toy Model

  28. 28.

Figure S13. Individual feature functions and target and SNN output for Toy
Model.

  29. 29.

Figure S14. Target and SNN output with noise for Toy Model.

##  S1 Text S1. SNN Training Method

Based on 15 single features, we modeled landslide susceptibility following the
SNN architecture. Our SNN model allowed for the exploration of multiple
combinations of parameters (e.g., 120 composite features for Level-2) without
relying on an expert’s choices, precondition, or classification of input
features. Our model found the most relevant features using a tournament
ranking approach we designed to quantify each control’s contribution to
susceptibility ($S_{j}$, where $j$ corresponds to single or composite
features) based on multiple steps of knowledge distillation. By superposing
$S_{j}$, we produced the total landslide susceptibility map, $S_{t}$, a number
ranging from 0 to 1 as the final product. Below are the details of SNN
training method.

###  S1.1 Algorithm

The flow diagram of the superposable neural networks (SNN) training method is
presented in Figure S1. The SNN is an additive model [Hooker (2004), S. Tan .
(2018)] with a unique architecture described by eq. (2) and Fig. 1, and a
unique training method explained here. Radial basis (Gaussian) activation
functions are used, where each neuron is connected to one input only. The
method consists of two main stages: feature-selection and knowledge
distillation. A set of optimal composite features are identified in the
feature selection stage, and multiple steps of knowledge distillation are used
to create an SNN. The SNN inputs are composite features corresponding to a
specific location, and the SNN output is the estimated landslide
susceptibility ($S_{t}$) at that location. A validation of our approach using
a toy model where we demonstrate the ability to extract boolean relationships
from synthetic data is included in Supporting Information.

The method can be summarized by the following steps:

  1. 1.

Multivariate polynomial expansion: composite features are generated.

  2. 2.

Tournament ranking: an automated feature selection method we have designed for
finding the features that are most relevant to the model.

  3. 3.

Multistage training (MST): a second-order deep learning technique for
generating a high-performance teacher network.

  4. 4.

Fractional knowledge distillation: a technique we designed for separating the
contribution of each feature to the final output.

  5. 5.

Parallel knowledge distillation: standard knowledge distillation individually
applied to networks corresponding to each feature.

  6. 6.

Network superposition: merging single layer networks corresponding to each
feature into one SNN.

These steps are explained in detail below.

###  S1.2 Multivariate polynomial expansion

Given $n$ features $x_{1},x_{2},..x_{n}$, we generate $M$ composite features
$\chi_{1},\chi_{2},..,\chi_{M}$ according to a predetermined maximum composite
feature level.

Ex. 1: If the original number of features is 3 and the maximum composite
feature level is Level-3, then we generate 13 composite features
$[\chi_{1},\chi_{2},..,\chi_{13}]=[x_{1},x_{2},x_{3},x_{1}*x_{2},x_{1}*x_{3},x_{2}*x_{3},x_{1}*x_{2}*x_{3},x_{1}^{2}*x_{2},x_{1}^{2}*x_{3},x_{2}^{2}*x_{1},x_{2}^{2}*x_{3},x_{3}^{2}*x_{1},x_{3}^{2}*x_{2}]$.

In this work, we have used 15 original features with a maximum composite
feature Level-2. Because Level-3 performs marginally better than Level-2, we
consider the Level-2 SNN as our optimal SNN. With 15 original features and the
maximum composite feature Level-2, we generate a total 120 composite features.
All features are standardized with zero-mean and unit-variance.

###  S1.3 Tournament ranking

Our feature selection technique is based on a point system and uses a
combination of backwards elimination and forward selection [L. Xu Zhang
(2001)] as building blocks. The composite features generated in the previous
steps are randomly arranged into groups, with each group containing a subset
of the features. Each feature group is used to train a simple neural network
model. After the network is trained, backwards elimination is applied using
area under the receiver operating characteristic curve (AUROC) as the
performance criterion (Supporting Information). The top performing feature in
the group receives a point. This process is repeated many times; several
thousand groups were generated in the training of each SNN in this work.
Features are ranked according to the points they accumulated. Forward
selection is then applied in the order of the feature ranking to select the
features that will be passed on to the next step.

The second-order Levenberg-Marquardt algorithm [Yu Wilamowski (2011)] was used
in training the individual neural networks models. It should be noted that
using second-order training is essential for the practicality of this step.
Unlike first-order training algorithms (based on gradient descent) that
require manual hyper parameter tuning, second-order training algorithms are
robust. In addition, second-order training can achieve better performance with
fewer parameters [Ozyildirim Kiran (2020), Le . (2011), Battiti (1992), H. Tan
Lim (2019), Montavon . (2012), Yu Wilamowski (2011), Wilamowski Yu (2010)].
This allows the automation of the process, and reduces the memory requirements
for training the networks, yielding a more efficient parallel implementation
on multicore processors.

###  S1.4 Multistage training

The high-ranked features that are passed on from the previous step are used to
train a high-performance DNN. We chose MST as our DNN model, since it has
shown superior performance in similar applications [Youssef . (2015), Bouchard
Youssef (2018), Youssef . (20181)].

###  S1.5 Fractional knowledge distillation

Knowledge distillation is a technique to reduce model complexity, by using the
soft output of a more complex teacher DNN as the target of a less complex
student DNN [Hinton . (2015)]. The MST in the previous step acts as our
teacher network.

We have designed a variation of knowledge distillation that allows us to
isolate the contribution of each feature to the estimated output. We call this
variation fractional knowledge distillation, a term that is inspired by the
fractional distillation technique in chemistry. We illustrate this using a
step-by-step example for the case of two features. This can be easily
generalized to any number of features.

Ex. 2: Assume that two composite features $[\chi_{1},\chi_{2}]$ are passed on
from the feature selection stage, and ordered according to importance where
$\chi_{1}$ is the most important. Let $ts_{0}$ be the set of soft targets
obtained from the MST output:

  1. 1.

Save a copy of $ts_{0}$, named $ts_{0c}$

  2. 2.

Train a simple DNN $net_{1,1}$ using only $\chi_{1}$ as input and $ts_{0}$ as
an output

  3. 3.

Obtain $o_{1,1}$, the set of outputs of $net_{1,1}$

  4. 4.

Update $ts_{0}$ to $ts_{0}-o_{1,1}$

  5. 5.

Train a simple DNN $net_{2,1}$ using only $\chi_{2}$ as input and $ts_{0}$ as
an output

  6. 6.

Obtain $o_{2,1}$, the set of outputs of $net_{2,1}$

  7. 7.

Update $ts_{0}$ to $ts_{0}-o_{2,1}$

  8. 8.

Evaluate performance by calculating AUROC using
$\sum_{i=1}^{2}\sum_{j=1}^{1}o_{i,j}$ and $ts_{0c}$

  9. 9.

Train a simple DNN $net_{1,2}$ using only $\chi_{1}$ as input and $ts_{0}$ as
an output

  10. 10.

Obtain $o_{1,2}$, the set of outputs of $net_{1,2}$

  11. 11.

Update $ts_{0}$ to $ts_{0}-o_{1,2}$

  12. 12.

Train a simple DNN $net_{2,2}$ using only $\chi_{2}$ as input and $ts_{0}$ as
an output

  13. 13.

Obtain $o_{2,2}$, the set of outputs of $net_{2,2}$

  14. 14.

Update $ts_{0}$ to $ts_{0}-o_{2,2}$

  15. 15.

Evaluate performance by calculating AUROC using
$\sum_{i=1}^{2}\sum_{j=2}^{2}o_{i,j}$ and $ts_{0c}$

  16. 16.

Repeat $n$ times until the performance stops improving

Each DNN above consists of only a few neurons and is trained for a small
number of epochs where the contribution of each feature is gradually
determined to avoid numerical instabilities. The number of neurons and epochs
are hyper parameters that can be tuned based on the data.

###  S1.6 Parallel knowledge distillation

The outputs from groups of networks, corresponding to each feature from the
previous step, are added together to yield one soft target per feature.
Knowledge distillation is separately used to train a single SNN layer for each
feature.

Ex. 3: Following the previous example:

  1. 1.

Create two soft targets: $ts_{1}=\sum_{j=1}^{n}o_{1,j}$, and
$ts_{2}=\sum_{j=1}^{n}o_{2,j}$

  2. 2.

Train a single layer network $net_{1}$ using $\chi_{1}$ as input and $ts_{1}$
as an output

  3. 3.

Train a single layer network $net_{2}$ using $\chi_{2}$ as input and $ts_{2}$
as an output

###  S1.7 Network superposition

The single layer networks from the previous step are merged together to create
the SNN, by adding an output layer that sums up the outputs of all the
networks from the previous step. The connection weights at the output layer
are set to one. The output of the SNN is a continuous value between 0 and 1,
which determines the network’s estimation of landslide susceptibility at a
specific location.

Ex. 4: Following the previous example, an SNN is created with $\chi_{1}$ and
$\chi_{2}$ as inputs and $O=o_{1}+o_{2}$ as the output, where $o_{1}$ is the
output of $net_{1}$ and $o_{2}$ is the output of $net_{2}$.

###  S1.8 Final notes

In this work, we have created three SNNs for three regions. The data samples
from each region were partitioned into roughly 70% for training and 30% for
testing. All reported results in the paper were obtained using the testing
portion of the data.

Class imbalance was taken into consideration when training the networks. Given
that the percentage of positive targets (locations containing a landslide) in
each region is significantly smaller than negative targets (locations with no
landslide), positive targets were weighted higher than negative targets in the
training cost functions following the approach in [Youssef . (20182)].

##  S2 Text S2. Description of Landslide Inventory

The landslide inventory was obtained from manual mapping supplemented by an
automatic detection technique. The basic procedure is as follows. We initially
mapped landslides using Landsat 8 satellite imagery from November 12, 2017. We
then used convolutional neural networks (CNN) to detect landslides
automatically, according to [Ghorbanzadeh . (2019), Prakash . (2020)]. The CNN
method used Landsat images and topographic and climatic attributes of mean
curvature, elevation, local relief, mean annual precipitation, slope, failure
index, and wetness to detect the landslides and test against the manually
mapped landslides. We then manually corrected landslides from the automatic
detection method using Landsat images, high-resolution satellite images from
Google Earth, and a 4-band Planetscope Scene with a 3 m resolution. Manual
correction is needed due to a common problem of inaccurately representing
landslide areas because of aggregated features and small detected features
that are not related to landslides in automatically mapped inventories [Marc
Hovius (2015), Parker . (2011)]. We divided aggregated features into multiple
landslides following suggestions from a previous study [Marc Hovius (2015)].
We used the manually corrected, automatically mapped landslides as our final
landslide inventory (hereafter called semi-automated landslides).

Most landslides are a combination of source and deposit, which are difficult
to distinguish in coarse resolution Landsat bands. Whenever possible, we
excluded debris transport or deposits and only mapped landslide scars
associated with source areas. In addition, we only assessed regions where
landslides have the potential to occur or be detectable. Thus, areas of
topographic slope less than 0.06 and alpine areas without vegetation cover
were excluded from our landslide mapping and analysis. A slope threshold of
0.06 was determined to be the minimum slope along which landslides occur based
on a cumulative distribution function of slope from observed landslides in the
eastern end of Himalaya. The criteria based on terrain characteristics such as
slope or local relief has been used in previous studies to constrain the area
of the landslide analysis [Parker . (2011)]. Alpine areas were classified
using spectral signatures representing snow cover in Landsat 8 imagery from
February 2018. The areas excluded from our mapping are shown in gray dashed
outlines in Fig. 2.

Our landslide inventory is comprised of semi-automated landslides includes
2,289 landslides covering a total area of 2.83$\times$107 m2, and a landslide
density of 0.007 over a 4.19$\times$109 m2 mapped region (Table S1). Areas
commissioned by semi-automatic detection but not manual mapping were
$\scriptstyle\sim$0.1, $\scriptstyle\sim$0.4, $\scriptstyle\sim$0.1% of the
total study area, while areas omitted by automated detection were
$\scriptstyle\sim$0.2, $\scriptstyle\sim$0.6, $\scriptstyle\sim$0.1% of the
total area.

The manually and semi-automatically detected landslides show a good
correspondence ($>$90% match for landslides $>$4 pixels (3,600 m2) based on
object identification that examine the existence of overlapping areas.
Generally, most landslides missing from the manually detected inventory are
objects with a small number of pixels that are not easily and objectively
detected by humans. Semi-automated landslides with $\leq$4 pixels comprise
$\sim$7.5% of total landslide areas. When comparing these pixels with high-
resolution Planetscope Scene satellite images during the post-processing
procedure, we find that many of these pixels are indeed small landslides
showing different spectral signatures. Thus, we include these semi-automatic
landslides with $\leq$ 4 pixels in our final inventory. The area frequency
distribution of our landslides from manual and semi-automatic mappings before
2017 shows a similar distribution to pre-2007 landslides from a nearby eastern
Himalaya region that are manually mapped using 15-30 m resolution ASTER and
Landsat images [Larsen . (2010), Larsen Montgomery (2012)]. Previous studies
[Larsen . (2010), Larsen Montgomery (2012)] used a landslide area of 100,000
m2 for the transition between soil and bedrock landslides in the nearby area.
Considering that $<$1% of landslides in numbers and $<$20% of total landslide
area are greater than 100,000 m2, most of the landslides in our semi-automated
inventory are likely soil landslides.

##  S3 Text S3. Feature Descriptions and Limitations

Here we explain how we obtained and calculated the 15 single features. The 15
single features include aspect (_Asp_), mean curvature (_Curv_ M), planform
curvature, profile curvature, total curvature, discharge, distance to channel
(_Dist_ C), distance to faults (_Dist_ F), distance to the Main Frontal Thrust
and suture zone (_Dist_ ${}_{\emph{MFT}}$), drainage area, elevation (_Elev_),
local relief (_Relief_), mean annual precipitation (_MAP_), number of extreme
rainfall events (_NEE_) and slope. The inclusion of these variables is based
on previous studies that examined the controls of landslides in the Himalayan
regions [Devkota . (2013), Regmi . (2014), Mandal Mandal (2018), Chowdhuri .
(2021)].

First, topographic variables such as slope, aspect, local relief, curvature,
distance to channel, and drainage area were calculated from a 90 m SRTM
Digital Elevation Model (DEM) [United States Geological Survey EarthExplorer
(2020)]. Although a higher-resolution 30 m DEM is available, it contains
missing values within our study area. Thus, we used a 90 m DEM for calculating
topographic variables. Slope was calculated as the steepest descent gradient
using an 8-direction (D8) flow routing method [Schwanghart Scherler (2014)].
We calculated aspect, the direction of slope face, as the angle in degrees
clockwise from north given by the components of the 3-D surface normal. The
surface normal was calculated using the $x$, $y$, and $z$ components of each
pixel. Local relief was calculated as the range in elevation within a 2.5 km
radius circular window. We used a 2.5 km radius window because it is similar
to the length scale of across-valley widths of channels in the range front
where most landslides are. Local relief at this scale allows us to quantify
the spatial variation of topographic relief relevant to landslides on these
fluvial valleys. Curvature was calculated as the second derivative of the 90 m
SRTM DEM. We calculated mean, planform, profile, and total curvatures using
TopoToolbox 2 [Schwanghart Scherler (2014), Schmidt . (2003)].

To calculate distance from channel, we first determined flow direction using
D8 flow routing. The flow direction was carved through topographic depressions
and flat areas to avoid sinks and generate a continuous drainage system. We
then imposed a minimum drainage area of 1 km2 needed to initiate stream before
extracting a stream network based on the flow direction. Using the stream
network, we calculated the distance of each pixel in the DEM to the nearest
location in the stream network.

We acquired _MAP_ and _NEE_ from a previous study [Bookhagen Burbank (2010)]
that analyzed the Tropical Rainfall Measuring Mission (TRMM) 2B31 datasets
from January 1998 to December 2009. Daily rainfall and _MAP_ values were
integrated from 3-hour measurements over 12 years. To calculate _NEE_ , the
90th percentile of daily rainfall total for each pixel was determined for the
12-year measurement period [Bookhagen Burbank (2010)]. Only days with measured
rainfall were included in calculating the probability density function. The
number of days per year with daily rainfall total above the 90th percentile
was counted as _NEE_ [Bookhagen Burbank (2010), Bookhagen (2010)]. The
resolution of the original _MAP_ and _NEE_ datasets in our study area was
$\sim$2.6 km, which we resampled to 30 m resolution to be consistent with the
resolution of our landslide inventory. To calculate the drainage area, we
first calculated D8 flow directions of stream networks and calculated the
number of upstream cells that contribute to each pixel. The number of cells
can then be converted into a drainage area. Discharge was calculated by
summing upstream contributing cells weighted by their _MAP_ to account for
spatially varying precipitation patterns. Using these weights, cells with
higher _MAP_ values will contribute more to total discharge than cells with
lower precipitation values.

Previous studies [Parker . (2011), C. Xu . (2014)] have shown that distance to
fault ruptures is a good predictor for the occurrence of earthquake-induced
landslides. We do not have information on active fault planes at depth and
ground peak acceleration patterns for past earthquakes in these regions. Thus,
we calculated _Dist_ ${}_{\emph{MFT}}$ for our study regions as each pixel’s
Euclidean distance from the closest point on traces of the Main Frontal Thrust
(MFT) and suture zones mapped by Taylor and Yin. [Taylor Yin (2009)] These
faults represent potentially active faults in our study area [Haproff .
(2019), Haproff . (2020)]. Because the suture zone is located far to the
north, _Dist_ ${}_{\emph{MFT}}$ largely reflects the distance to the MFT. In
addition, we calculated _Dist_ F as the Euclidean distance of each pixel from
boundaries separating all lithologic units reported in [Haproff . (2019)]. We
include _Dist_ F because bedrock tends to be more damaged near major
lithologic boundaries due to faulting, which may influence landslide
occurrences. The Euclidean distance was calculated using ArcGIS 10.6.

We note that there are limitations of our data from the eastern end of
Himalaya. Landslide and input feature data have relatively coarse spatial
resolutions and are based on limited spatial and temporal information (e.g.,
30 m resolution Landsat satellite images from 2017 [United States Geological
Survey EarthExplorer (2020)], 90 m resolution SRTM DEM [United States
Geological Survey EarthExplorer (2020)], and $\scriptstyle\sim$2.6 km
resolution TRMM data over 12 years [Bookhagen Burbank (2010)]). We don’t have
readily available high-resolution data on topography, surface materials (e.g.,
soil depth, bedrock structures, lithology), and climatic and ecohydrologic
conditions (e.g., landslide-triggering storm intensity, time-series
precipitation intensity, vegetation types). Though limited high-resolution
satellite imagery is available (e.g., Planetscope Scene) in the area, we used
relatively coarse 30 m resolution Landsat images because 1) Landsat images are
globally available, open-source satellite images with a
$\scriptstyle\sim$40-year historic archive, 2) topographic, climatic, and
geologic feature data have coarser resolutions than 30 m, and 3) we cover a
large region of the eastern Himalaya (a total area of 4.19$\times$109 m2,
4.66$\times$106 pixels at 30 m). For modeling a regional-scale study covering
a large area, limited input data resolution, and computational costs, the use
of 30 m resolution for our model was inevitable. Although our inventory is
based on coarse 30 m resolution Landsat images, our inventory generated by
semi-automatic mapping includes many small landslides missed by humans, but
detected by an automatic detection algorithm based on a DNN and confirmed by
high-resolution satellite images (see section S1). Thus, our landslide
inventory captures the regional-scale, spatial distributions of landslide
occurrences and provides essential information for regional-scale landslide
susceptibility models. Due to the uncertainties in feature data, our landslide
susceptibility model and feature analysis may also contain similar
uncertainties and may change with updated, high resolution input data.
However, our work presents more accurate and highly interpretable landslide
susceptibility model based on currently available data, which can be readily
available for hazard prediction and mitigation. Our SNN will be also easily
applicable to other locations and different datasets when global, open-source,
high-resolution datasets and high-performance computational resources become
more available in the future.

##  S4 Text S4. Explanations for Landslide Controls

Here we explain how key results from the SNN can provide insights into the
physical controls and mechanisms of landslides. A physically-derived metric,
failure index (_FI_), based on limit equilibrium represents landslide
propensity on an infinite slope for cohesionless material considering the pore
pressure effect [Montgomery Dietrich (1994), Moon . (2011)]. This metric is
similar to an inverse of factor of safety and can be related to total
landslide susceptibility ($S_{t}$). We can formulate _FI_ as:

$\displaystyle
FI=\frac{S}{S_{0}}\left(1-W\frac{\rho_{w}}{\rho_{s}}\right)^{-1}$ (2)

where $S_{0}$ is the threshold slope, $S$ is the local slope, $\rho_{s}$ is
the wet bulk density of soil (2.0 g/cm3), $\rho_{w}$ is the bulk density of
water (1.0 g/cm3), and $W$ is wetness that represents the degree of subsurface
saturation. $W$ varies from 0 (unsaturated) to 1 (fully saturated). This
equation can be further expanded by Taylor series:

$\displaystyle
FI\approx\frac{S}{S_{0}}\left(1+W\frac{\rho_{w}}{\rho_{s}}+W^{2}(\frac{\rho_{w}}{\rho_{s}})^{2}+O(W^{3})\right)=\frac{S}{S_{0}}\left(k(W)\right).$
(3)

The approximated _FI_ has three components, local slope $S$, threshold slope
$S_{0}$, and _k(W)_ that represents the degrees at which landslides are
promoted by subsurface saturation. _k(W)_ varies from 1 (unsaturated) to 2
(fully saturated).

A change in climatic conditions can increase volumetric water content and
porewater pressure and lead to an increased degree of subsurface saturation
(i.e., $W$) and subsequently induce slope failure. Previous physically-based
slope stability models consider various climatic factors (e.g., rainfall
amount and intensity, subsurface convergence flow) to deduce the degree of
subsurface saturation to model rainfall-induced landslide occurrences
[Montgomery Dietrich (1994), Baum . (2002), Baum . (2010)]. For example, the
Shallow Landslide Stability model (SHALSTAB) [Montgomery Dietrich (1994),
Dietrich . (2001)] uses the topographic wetness index, proposed by Beven and
Kirkby (1979) [Beven Kirkby (1979)], to calculate subsurface saturation
considering the convergence of shallow subsurface flow from up-slope drainage
areas for a given steady-state precipitation. On the other hand, the Transient
Rainfall Infiltration and Grid based Regional Slope stability model (TRIGRS)
[Baum . (2002), Baum . (2010)] calculates transient pore pressure development
due to vertical rainfall infiltration from rainfall intensity.

In reality, both subsurface convergence and rainfall infiltration are
essential contributors to subsurface saturation and need to be implemented in
physically-based slope stability models. However, measuring precipitation
intensity, moisture availability, or subsurface convergence and saturation is
difficult, especially in rural mountainous areas with limited accessibility.
Thus, the identification and quantification of various climatic contributions
that promote landslides is challenging.

According to the Level-2 SNN model results, the most important,
differentiating features for landslides in this area are the product of slope
and mean annual precipitation (_MAP_) (N-S region) or the product of slope and
number of extreme events (_NEE_) (NW-SE and E-W regions). This result implies
that local precipitation infiltration influenced by precipitation rate and
intensity, represented by _MAP_ and _NEE_ , may serve as a first-order control
on _k(W)_ in equation 2. The absence of drainage area or discharge as a
contributing feature to susceptibility may suggest that subsurface flow
convergence may be a second order contributor to the landslides in the eastern
end of Himalaya. However, we cannot rule out the possibility that the
importance of topographic convergence was masked due to the low-resolution of
input topographic and rainfall data, e.g. [Leonarduzzi . (2021)]. These
factors can be further examined in future studies using high-resolution
topographic and climate data.

In addition, the SNN identified aspect, the direction of slope face, as
another primary feature that influences landslide occurrences. Previous
studies considered slope aspect preference in terms of: 1) vegetation activity
that affects root cohesion [McGuire . (2016)], or 2) the orientation of wind-
driven rainfall. To examine vegetation activity across hillslope aspect, we
calculated the normalized difference vegetation index (_NDVI_) following the
USGS procedure [Zanter (2016)]. We first converted Landsat 8 Level-1 Digital
Numbers to top-of-atmosphere (TOA) reflectance. TOA reflectance eliminates the
impact of different solar angles and illumination geometries and is calculated
as:

$\displaystyle\rho\lambda=\frac{M_{\lambda}Q_{cal}+A{\rho}}{\cos(\theta_{SZ})},$
(4)

where $\rho\lambda$ is the TOA reflectance, $M_{\lambda}$ is the band-specific
multiplicative rescaling factor from the Landsat 8 metadata, _Q_ cal is the
standard product pixel values, $A\rho$ is the band-specific additive rescaling
factor from the metadata, and $\theta_{SZ}$ is the local solar zenith angle.

We use the corrected bands 4 and 5 from Landsat 8 to calculate NDVI as:

$\displaystyle\emph{NDVI}=\frac{\mbox{band}\,5-\mbox{band}\,4}{\mbox{band}\,5+\mbox{band}\,4},$
(5)

where bands 4 and 5 represent the visible and near-infrared light reflected by
vegetation, respectively. Healthy vegetation with high photosynthetic capacity
absorbs a larger proportion of incident visible light while reflecting a
greater portion of near-infrared light compared to sparse or unhealthy
vegetation [Tucker . (1986)]. Therefore, an _NDVI_ value close to 1 suggests a
higher density of healthy vegetation and green leaves while a value near 0
might indicate unhealthy or no vegetation. We utilized Landsat 8 satellite
imagery from October 2015, November 2017, and February 2018 [United States
Geological Survey EarthExplorer (2020)] for our analyses of _NDVI_. These
months were selected to characterize _NDVI_ values before and after the summer
monsoon season, during which a large proportion of landslides are suspected to
occur because of intense rainfall. We excluded summer months from our analyses
because of the abundant cloud cover present in those images, which masks the
visibility of the land surface.

_NDVI_ plotted against aspect in our study area shows a broad distribution of
high values centered around values corresponding to south-facing slopes.
However, this _NDVI_ distribution is different from the observed peak of _S
Asp_ around 145°to 180°(Figure S11). This result may imply that more
landslides on south-facing slopes are likely due to orographic precipitation
patterns caused by moisture delivery from the south rather than through the
effects of vegetation. Previous work has characterized the northward moisture
transfer to this study area from the Bay of Bengal during monsoon seasons
[Bookhagen Burbank (2010), Barros . (2004), Yang . (2018)]. Thus, we think
that the SNN-identified primary feature of aspect supports the influence of
aspect-related differences in microclimate (e.g., moisture availability) on
landslide occurrences in this area.

##  S5 Text S5. Assessments of Model Performances

We evaluate the performance of the SNN compared to traditional approaches
using several performance metrics, including the area under the receiver
operating characteristic curve (AUROC), the area under the success rate curve
(AUC of success rate), accuracy, sensitivity (i.e., probability of detection,
POD), specificity (i.e., probability of false detection, POFD), and POD-POFD
following the literature. See, for example [Prakash . (2020)]. AUROC and AUC
of success rate are cutoff-independent performance criteria while accuracy,
POD, and POFD are cutoff-dependent. The AUROC is calculated as the area under
a curve created by plotting the true positive rate against the false positive
rate at various thresholds along a feature’s range. The AUC of success rate is
calculated by plotting the fraction of landslide area correctly identified
against the fraction of total area classified as positive (i.e., predicted
landslides). AUC ranges between 0 and 1, with 1 indicating a perfect
classifier and 0.5 indicating a random model. After generating a threshold-
modeled landslide map based within the 30% testing partition using the optimal
threshold corresponding to the point closest to [0,1] on an AUROC curve,
accuracy is calculated as the fraction of landslide and non-landslide area
correctly classified by the model relative to all studied areas. POD and POFD
measure the proportion of landslide areas correctly classified relative to all
observed landslide areas and the proportion of incorrectly classified
landslide areas relative to all observed non-landslide areas, respectively.

We calculate these metrics for all 15 single original features, a physically-
based slope stability model (SHALSTAB), two statistical methods (logistic
regression and likelihood ratios), a first-order artificial neural network
(ANN) based on the adaptive moment estimation (Adam) optimization algorithm,
and Level-1 and Level-2 SNNs.

First, we investigated each of the 15 single features as individual
classifiers for landslide occurrences (Table S3). Second, we assessed the
propensity of landslides using a topographic metric called the failure index.
The failure index (_FI_) is the ratio of driving to resisting force on a
hillslope, which is the inverse of the factor-of-safety. _FI_ is modified from
the SHALSTAB model that couples infinite slope stability and steady-state
hydrology for cohesionless material [Montgomery Dietrich (1994), Dietrich .
(1995), Dietrich . (2001), Moon . (2011)]. Considering that landslides smaller
than 100,000 m2 (the upper bound for soil landslides in this region [Larsen .
(2010), Larsen Montgomery (2012)]) constitute $>$99% of landslides in number
and $\scriptstyle\sim$80% of total landslide area, we assume that most
landslides in our inventory are soil landslides.

To calculate the _FI_ , we first determine the spatial distribution of wetness
($W$), which represents the degree of subsurface saturation. $W$ is calculated
as a ratio between local hydraulic flux from a given steady-state
precipitation relative to that of soil profile saturation [Montgomery Dietrich
(1994)]:

$\displaystyle W=\frac{h}{z}=\frac{qA}{bT\sin\theta}$ (6)

where $h$ is the saturated height of the soil column ($L$), $z$ is the total
height of the soil column ($L$), $q$ is the steady-state precipitation during
a storm event ($L$/$T$), $A$ is the drainage area (_L 2_) draining across the
contour length $b$ ($L$), $T$ is the soil transmissivity when saturated
($L\textsuperscript{2}/T$), and $\theta$ is the local slope in degrees. $W$
varies from 0 (unsaturated) to a capped value of 1 (fully saturated). We used
a base value of 1$\times$10-4 m2/s for $T$ following Moon et al. [Moon .
(2011)] We used the spatial distribution of _MAP_ [Bookhagen Burbank (2010)]
to represent the steady-state precipitation, $q$. We then calculated the
spatial distribution of _FI_ as:

$\displaystyle
FI=\frac{S}{S_{0}}\left(1-W\frac{\rho_{w}}{\rho_{s}}\right)^{-1}$ (7)

where $S_{0}$ is the threshold slope set at 45°, $S$ is the local slope,
$\rho_{s}$ is the wet bulk density of soil (2.0 g/cm3), and $\rho_{w}$ is the
bulk density of water (1.0 g/cm3). To examine whether the model performance of
_FI_ is different when predicting all landslides vs soil landslides, we
include the results of the performance metrics for _FI_ calculated using all
landslides and soil landslides in Table S3.

Third, we applied two statistical models, logistic regression and likelihood
ratios, to assess landslide susceptibility. Logistic regression (hereafter,
LogR) is based on a multivariate regression between a binary response of
landslide occurrence and a set of predicting features that are continuous,
discrete, or a combination of both types [Lee (2005)]. To build these models,
we consider only one curvature metric following Lee [Lee (2005)], instead of
using all four different curvatures. We selected _Curv M_ to build the
statistical models. In addition, we consider log10(drainage area) and
log10(discharge) because of their inverse power-law relationships with
landslide and debris flow incision [Stock Dietrich (2003), Stock Dietrich
(2006)]. The relationship between features and landslide occurrence can be
displayed as:

$\displaystyle p=\frac{e^{c}}{e^{c}+1}$ (8)

where _p_ is the probability of landslide occurrence that varies from 0 to 1
in an S-shaped curve, and _c_ is the linear combination of features:

$\displaystyle c=b_{0}+b_{1}x_{1}+b_{2}x_{2}+...+b_{n}x_{n}$ (9)

where $x_{i}$ ($i$ = 1, 2, . . ., $n$) represents each feature, $b_{i}$
represents the optimized coefficient, and $b_{0}$ represents the intercept of
the model. Utilizing Eqs. 8 and 12, we obtain an extended expression for the
LogR model relating the probability of landslide occurrence $p$ and multiple
features:

$\displaystyle\mbox{logit}(p)=\log\left(\frac{p}{1-p}\right)=b_{0}+b_{1}x_{1}+b_{2}x_{2}+...+b_{n}x_{n}$
(10)

where $\log$ is the natural log. To determine any possible collinearity
between features, we calculate the correlation coefficient ($R$) between all
combinations of 12 features (Table S6). We observed maximum absolute values of
$R$ = -0.828 (N-S), 0.717 (NW-SE), and 0.857 (E-W), which are below the
threshold of 0.894 corresponding to a variance inflation factor of <5,
indicating low collinearity between features [Stine (1995), Kavzoglu . (2014)]
and thus we used all 12 features. We treated aspect as a discrete feature due
to its nonlinear relation with landslide occurrences. The best-fit coefficient
values are shown in Table S7.

Similar to the SNN, the LogR method provides information on the importance of
variables through the best-fit coefficients. To compare those results, we
determined top features that differentiate areas with and without threshold
modeled landslides based on the LogR output following similar procedures that
we used for the SNN. The output of LogR ranges between 1.06$\times$10-6 to
0.820 for the N-S region. The threshold value ($t$) of 0.005 that corresponds
to the point closest to [0,1] on an ROC curve (i.e., a perfect classifier) was
used to classify landslide (_ls_) and non-landslide areas (_nls_) for the N-S
region. We calculated $\Delta\bar{cn}_{i}$ as the difference between the
average value of a feature multiplied by its respective coefficient for _ls_
and _nls_ areas, then divided by an adjusted threshold that was transformed
from $t$ (i.e., 0.005) from the LogR output according to the equation below:

$\displaystyle\Delta\bar{cn}_{i}=\frac{\bar{c}_{i,ls}-\bar{c}_{i,nls}}{t_{a}}\,\mbox{~{}where~{}}\,t_{a}=\log\left(\frac{t}{1-t}\right)-b_{0},$
(11)

where $\bar{c}_{i,a}$ is the average feature $i$ value multiplied by its
respective coefficient for areas ($a$) of _ls_ or _nls_ , and $t_{a}$ is the
adjusted threshold value based on $b_{0}$, the overall intercept value
determined by the LogR model, and $t$, the threshold determined using the ROC
curve (i.e., 0.005 for the N-S region). We transformed $t$ to $t_{a}$ and used
it for normalization to enable the direct comparison of results between LogR
and SNN. For the SNN-determined primary features, we calculate
$\Delta\bar{Sn}_{j}$ as $\Delta\bar{S}_{j}$ divided by the threshold that
classify landslides (i.e., 0.406 for the N-S regions). The value of 1 in both
$\bar{Sn}$ and $\bar{cn}$ represent the threshold susceptibility that
classifies _ls_ and _nls_ areas. The results of the identified primary
controls of landslides that induce large differences in average susceptibility
between landslide and non-landslide areas in the N-S regions from LogR, SNN
Level-1, and SNN Level-2 are shown in Figure S7. All methods identified
climate-related factors (e.g., MAP, NEE, Asp) as primary controls; however,
only SNN Level-2 was able to identify the importance of the composite feature
_MAP_ *_Slope_.

The likelihood ratio method uses the relationship between observed landslide
occurrences and controlling feature ranges. Previous studies have quantified
the ratio of the probability of landslide occurrences to the probability of
non-occurrences or all-occurrences and referred to it as the likelihood ratio,
frequency ratio, or probability ratio [Lee (2005), Akgun (2012), Reichenbach .
(2018)]. In this study, we calculate likelihood ratios ($LR$) as the ratio of
the percentage of landslide pixels relative to total landslide pixels divided
by the percentage of pixels relative to the total area for the specific range
of feature values [Lee (2005), Akgun (2012)]. Landslide susceptibility for
each pixel is calculated as the sum of the corresponding $LR$ from each
feature’s value. A ratio of 1 and $>$1 indicates the average and above-average
likelihood of landslide occurrence within the feature range compared to that
of the study area. Conversely, values less than 1 indicate a below-average
likelihood. In this study, we used all 15 features with each feature’s range
divided into ten bins to calculate $LR$ and landslide susceptibility. The
first and last bins represent areas less than and greater than the 10th and
90th percentile of $LR$, respectively, with values between these bins split
into eight equal bin ranges.

Lastly, we applied a landslide susceptibility model based on an ANN using the
Adam optimizer. We prepared landslide feature datasets following [Tien Bui .
(2012)]. All features, except for planform, profile, and total curvatures,
used to build the SNN were used to apply the Adam optimizer. We divided each
of the 12 feature maps into 10 categories representing intervals for every
10th percentile. Then, we gave each category a normalized category value using
the following Max-Min normalization formula: [Fernandes Lona (2005)]

$\displaystyle v^{\prime}=\frac{v-\min(v)}{\max(v)-\min(v)}(U-L)+L$ (12)

where $v’$ is the normalized category value, $v$ is the attribute value, and
$U$ and $L$ are the upper and lower normalization range boundaries set as 0.1
to 0.9. Landslide inventories and feature maps were then partitioned into
training (70%) and testing (30%) regions using the same checkerboard pattern
applied when building the SNN.

The ANN analysis utilized 12 feature inputs and one network output. The
structure was composed of an input layer, three fully connected layers with
200, 200, and 2 neurons per layer respectively, a softmax layer, and a
classification output layer. The mini-batch size used for each training
iteration was set to 64. The maximum number of epochs was set to 200. The
initial learning rate was set to 2$\times$10-4. We used a piecewise method for
learning rate with the default learning rate drop factor of 0.1 that is
applied to the global learning rate every five epochs passed. The training
data was shuffled before every training epoch. An early stoppage of 10 epochs
without improvement in validation error was used to prevent over-fitting.

We determined a 95% confidence intervals of mean AUROC by conducting a 10-fold
cross validation for all statistical and neural network models utilized in
this study. We randomly and uniformly selected 50% of the original training
pixels from the checkerboard pattern partition and trained the model on this
portion of the training dataset. Subsequently, we tested the trained model on
50% of the testing dataset that was also selected randomly and uniformly. We
then calculated the AUROC for each trial. This procedure was then repeated 10
times for each method and the results were used to calculate the 95%
confidence interval for the mean AUROC and 2$\sigma$ range of AUROC from 10
validation tests (Table S5).

Our model assessments for single features indicate that _MAP_ [AUROC = 0.756
(N-S region)] and slope [AUROC = 0.696 (NW-SE), 0.760 (E-W)] are the highest
performing single features. The SNN produces a 19-26% improvement in AUROC
compared to a physically-based landslide model (e.g., failure index for all
landslides or soil landslides). Additionally, the SNN produced up to
$\scriptstyle\sim$8%, $\scriptstyle\sim$10%, and $\scriptstyle\sim$12%
increases in performance compared to the logistic regression, likelihood ratio
methods, and ANN-based susceptibility model, respectively. The physically-
based model of _FI_ produces slightly different AUROC for predicting all
landslides vs. soil landslides, but both AUROC values were lower than that of
the SNN (Table S3). Further investigation using performance metrics including
the AUC of the ROC, success rate, accuracy, POD, POFD, and POD-POFD reveals
that the SNN outperformed the tested ANN, and statistical and physical models
across all metrics (Table S4).

Although both likelihood ratios and the SNN reveal functional relationships
between susceptibility and features, the SNN is more useful for landslide
susceptibility assessment because it decouples individual feature
contributions and quantifies absolute contributions from features and feature
interdependencies. The SNN reveals the importance of features coupling slope
and climatic factors (e.g., _MAP*slope_ , _NEE*slope_) as a primary driver for
landslide occurrence. Accounting for these under-appreciated features and
feature interdependencies that are not generally implemented in statistical
methods or physically-based models can lead to a substantial increase in
performance.

##  S6 Text S6. Toy Model

There can be many solutions of models that can fit a dataset generated by
another model with varying degrees of accuracy. In order to test whether our
method can find a solution that closely matches a known existing model,
consider the following constrained toy application by generating a dataset
that represents a logical relationship and testing the behavior of our
algorithm.

Take the equation

$y=x1*x2+x3*x4-2*x1*x2*x3*x4,$ (13)

where $x1,x2,x3$ and $x4$ are Boolean values. It is easy to check that the
equation represents the logical relationship

$y=(x1\wedge x2)\vee(x3\wedge x4)\wedge(\overline{x1\wedge x2\wedge x3\wedge
x4})$

where the truth table is shown in Table S8.

We generated 1,000 random realizations of $x1_{n},x2_{n},x3_{n}$ and $x4_{n}$
and we calculate the corresponding value $y_{n}$ for each of these
realizations, where $n=1:1000$. We tested our algorithm by training an SNN
using $x1_{n},x2_{n},x3_{n}$ and $x4_{n}$ as the input and $y_{n}$ as the
target output, so as to test whether our method can infer the logical
relationship from the basic components using only the generated data samples.
Up to Level-4 composite features were used in this analysis, for a total
number of 15 features as shown in Table S9.

The resulting feature ranking (Figure S12) shows that our algorithm was able
to successfully isolate the composite features that exist in the relationship.
The truth table (Table S8) reveals that the higher ranking given to
$\\{x1*x2*x3*x4\\}$ corresponds to the fact that this feature can decisively
explain 50% of the logical relationship independently from the other features.
If $\\{x1*x2*x3*x4\\}$ is one, then $y$ is always zero. On the other hand, the
other features cannot decisively determine any part of the outcome on their
own, but they can decisively determine the outcome if they depend on
$\\{x1*x2*x3*x4\\}$.

Our method was able to find a solution to the logical relationship that
accurately matches the target output (Figure S13). Our model is given by

$\tilde{y}=f(x1*x2)+f(x3*x4)+f(x1*x2*x3*x4),$

where

$f(x1*x2)\approx\begin{cases}0.1,&\mbox{for~{}}x_{1}*x_{2}=0\\\
1.1,&\mbox{for~{}}x_{1}*x_{2}=1\end{cases}$

$f(x3*x4)\approx\begin{cases}-0.25,&\mbox{for~{}}x_{3}*x_{4}=0\\\
0.75,&\mbox{for~{}}x_{3}*x_{4}=1\end{cases}$

$f(x1*x2*x3*x4)\approx\begin{cases}0.15,&\mbox{for~{}}x_{1}*x_{2}*x_{3}*x_{4}=0\\\
-1.85,&\mbox{for~{}}x_{1}*x_{2}*x_{3}*x_{4}=1\end{cases}$

Note that subtracting the two ends of the function of each feature returns the
coefficient values in the original equation (Eq. 13) for each corresponding
variable:

$f(x1*x2):1.1-0.1=1,$ $f(x3*x4):0.75-(-0.25)=1,$
$f(x1*x2*x3*x4):-1.85-0.15=-2.$

Furthermore, the behavior of our method proved to be robust to noise. Figure
S14 demonstrates the results of the same experiment, only this time the data
was contaminated by adding Gaussian noise to $x1_{n},x2_{n},x3_{n}$ and
$x4_{n}$ prior to training and testing. Although the SNN output became
noisier, the noise effect can easily be removed with thresholding and does not
affect the outcome.

Fig. S1:  Superposable neural network training flow diagram. The details of
methodology are explained in the Data & Methods section (main text).

Fig. S2:  Elevation map of the eastern end of Himalaya. Yellow boxes indicate
our three N-S, NW-SE, and E-W oriented study regions. The inset map shows the
Himalayan region with our study area shown in a yellow box and national
borders shown with dark gray lines.

Fig. S3:  Spatial distribution of 15 features used in the superposable neural
network model. The 15 single features include (a) aspect (_Asp_), (b) mean
curvature (_Curv_ M), (c) planform curvature (_Curv_ Plan), (d) profile
curvature (_Curv_ Prof), (e) total curvature (_Curv_ Tot), (f) discharge, (g)
drainage area, (h) distance to channel (_Dist_ C), (i) distance to faults
(_Dist_ F), (j) distance to the Main Frontal Thrust and suture zone (_Dist_
${}_{\emph{MFT}}$), (k) elevation (_Elev_), (l) local relief (_Relief_), (m)
mean annual precipitation (_MAP_), (n) number of extreme events (_NEE_), and
(o) slope. Dashed lines mark the overlapping area between the N-S and NW-SE
region. Features in (e, f, g) are displayed on logarithmic scales.

Fig. S4:  Comparison between the feature ranges of distance to main faults
(_Dist MFT_), distance to faults (_Dist F_), number of extreme rainfall events
(_NEE_), mean annual precipitation (_MAP_), and aspect (_Asp_), each
normalized by the maximum feature value across all three regions. Red center
lines represent the median and top and bottom ends of the box represent the
25th and 75th quartiles, respectively. The ends of the dashed lines extending
from each side of the box plot represent 1.5 times the interquartile range or
the minimum or maximum values. Outliers are not shown in this figure. On the
$x$-axis, 1, 2, and 3 correspond to the N-S, NW-SE, and E-W regions,
respectively.

Fig. S5:  Spatial distribution of landslides, susceptibility, and important
controls for the NW-SE region. Maps are shown for (a) manually mapped
landslides in red polygons, (b) total susceptibility ($S_{t}$) from the
Level-2 superposable neural network, and (c) primary features identified as
locally important controls of landslides. Bar charts in (d) represent
$\Delta\bar{S}_{j}$ in descending order, and pie charts in (e) represent
average $S_{j}$ ($\bar{S}_{j}$) contributions to $ld$ and $nld$ areas for the
entire region. Details are the same as in Figure 2.

Fig. S6:  Spatial distribution of landslides, susceptibility, and important
controls for the E-W region. Maps are shown for (a) manually mapped landslides
in red polygons, (b) total susceptibility ($S_{t}$) from the Level-2
superposable neural network, and (c) primary features identified as locally
important controls of landslides. Bar charts in (d) represent
$\Delta\bar{S}_{j}$ in descending order, and pie charts in (e) represent
average $S_{j}$ ($\bar{S}_{j}$) contributions to $ld$ and $nld$ areas for the
entire region. Details are the same as in Figure 2.

Fig. S7:  Bar charts representing $\Delta\bar{Sn}_{j}$ for the (a) SNN - Level
2 and (b) SNN- Level 1 and $\Delta\bar{cn}_{i}$ for (c) logistic regression,
arranged in descending order. Details on the calculations of
$\Delta\bar{Sn}_{j}$ and $\Delta\bar{cn}_{i}$ are provided in the Supporting
Information, page S5. Features related to topography, aspect, climate, and
geology are shown in green, pink, blue, and brown or combinations thereof,
respectively.

Fig. S8:  Local analysis of climate vs slope susceptibility contributions to
landslide occurrence for the (a) N-S, (b) NW-SE, and (c) E-W study regions.
Values greater than 1 indicate a larger contribution from climate features
($\Sigma\Delta\bar{S}\textsubscript{\emph{climate}}$, $j$ = _Asp, NEE, MAP_)
relative to that of slope ($\Delta\bar{S}\textsubscript{\emph{Slope}}$).
Windows of no data are due to no modeled or mapped landslides.

Fig. S9:  Primary landslide controls of a) number of extreme events (_NEE_) *
slope, and b) slope * aspect (_Asp_) for the NW-SE region are shown. Details
are the same as in Figure 3.

Fig. S10:  Primary landslide controls of a) number of extreme events (_NEE_) * slope, b) aspect (_Asp_) * _NEE_ , and c) mean annual precipitation (_MAP_) * distance to main faults (_Dist MFT_) for the E-W region are shown. Details are the same as in Figure 3. Fig. S11:  The relationship between aspect and normalized difference vegetation index (_NDVI_), shown as thin lines, and $S\textsubscript{\emph{Asp}}$ from SNN level-1, shown as thick lines. Colors correspond to different regions while symbols shown as thin lines correspond to different times of measurements (October 2015, November 2017, and February 2018). Symbols on thin lines represent the averaged _NDVI_ value for a 20°, interval of aspect. Fig. S12:  Feature ranking. Fig. S13:  (A) individual feature functions. (B) Target output v.s. SNN output. Fig. S14:  Target output vs SNN output when data is contaminated with noise. Table S1:  Description of Landslide Inventory. Table S2:  Description and Ranges of 15 Features. Table S3:  AUROC of Models and Single Features. Entries in italics have negative correlations with landslide propensity AUROC of physical models and single features were normalized to be from 0 to 1 Based on 30% testing partition. Table S4:  Artificial Neural Network and Statistical Model Confidence Intervals. Table S5:  Performance Metrics. Table S6:  Correlation Metrics Between Features. Table S7:  Correlation Metrics Between Features. $x1$ | $x2$ | $x3$ | $x4$ | $x1*x2$ | $x3*x4$ | $x1*x2*x3*x4$ | $y$
---|---|---|---|---|---|---|---
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
0 | 0 | 0 | 1 | 0 | 0 | 0 | 0
0 | 0 | 1 | 0 | 0 | 0 | 0 | 0
0 | 0 | 1 | 1 | 0 | 1 | 0 | 1
0 | 1 | 0 | 0 | 0 | 0 | 0 | 0
0 | 1 | 0 | 1 | 0 | 0 | 0 | 0
0 | 1 | 1 | 0 | 0 | 0 | 0 | 0
0 | 1 | 1 | 1 | 0 | 1 | 0 | 1
1 | 0 | 0 | 0 | 0 | 0 | 0 | 0
1 | 0 | 0 | 1 | 0 | 0 | 0 | 0
1 | 0 | 1 | 0 | 0 | 0 | 0 | 0
1 | 0 | 1 | 1 | 0 | 1 | 0 | 1
1 | 1 | 0 | 0 | 1 | 0 | 0 | 1
1 | 1 | 0 | 1 | 1 | 0 | 0 | 1
1 | 1 | 1 | 0 | 1 | 0 | 0 | 1
1 | 1 | 1 | 1 | 1 | 1 | 1 | 0
Table S8:  Truth table. Feature | Level
---|---
$x1$ | 1
$x2$ | 1
$x3$ | 1
$x4$ | 1
$x1*x2$ | 2
$x1*x3$ | 2
$x1*x4$ | 2
$x2*x3$ | 2
$x2*x4$ | 2
$x3*x4$ | 2
$x1*x2*x3$ | 3
$x1*x2*x4$ | 3
$x1*x3*x4$ | 3
$x2*x3*x4$ | 3
$x1*x2*x3*x4$ | 4
Table S9:  Composite features.

## References

  * Adadi Berrada (2018) refx2Adadi, A. Berrada, M.  2018\.  Peeking inside the black-box: A survey on Explainable Artificial Intelligence (XAI) Peeking inside the black-box: A survey on explainable artificial intelligence (xai).  IEEE Access652138–52160.
  * Agarwal . (2020) agarwal2020neuralAgarwal, R., Frosst, N., Zhang, X., Caruana, R. Hinton, G.  2020\.  Neural additive models: Interpretable machine learning with neural nets Neural additive models: Interpretable machine learning with neural nets.  arXiv preprint arXiv:2004.13912.
  * Akgun (2012) akgun2012comparisonAkgun, A.  2012\.  A comparison of landslide susceptibility maps produced by logistic regression, multi-criteria decision, and likelihood ratio methods: a case study at İzmir, Turkey A comparison of landslide susceptibility maps produced by logistic regression, multi-criteria decision, and likelihood ratio methods: a case study at i̇zmir, turkey.  Landslides9193–106.
  * Barros . (2004) Barros2004Barros, A., Kim, G., Williams, E. Nesbitt, S.  2004\.  Probing orographic controls in the Himalayas during the monsoon using satellite imagery Probing orographic controls in the Himalayas during the monsoon using satellite imagery.  Natural Hazards and Earth System Science, Copernicus Publications on behalf of the European Geosciences Union4.
  * Battiti (1992) khaled_r3Battiti, R.  1992\.  First- and Second-Order Methods for Learning: Between Steepest Descent and Newton’s Method First- and second-order methods for learning: Between steepest descent and newton’s method.  Neural Computation42141-166.  10.1162/neco.1992.4.2.141
  * Baum . (2010) Baum2010Baum, R., Godt, J. Savage, W.  2010\.  Estimating the timing and location of shallow rainfall‐induced landslides using a model for transient, unsaturated infiltration Estimating the timing and location of shallow rainfall‐induced landslides using a model for transient, unsaturated infiltration.  Journal of Geophysical Research: Earth Surface115F3.
  * Baum . (2002) Baum2002Baum, R., Savage, W. Godt, J.  2002\.  TRIGRS—a Fortran program for transient rainfall infiltration and grid-based regional slope-stability analysis TRIGRS—a Fortran program for transient rainfall infiltration and grid-based regional slope-stability analysis.  US geological survey open-file report42438.
  * Ben-Menahem . (1974) Ben-Menahem1974Ben-Menahem, A., Aboodi, E. Schild, R.  1974\.  The source of the great Assam earthquake—an interplate wedge motion The source of the great Assam earthquake—an interplate wedge motion.  Physics of the Earth and Planetary Interiors94265-289.
  * Bergen . (2019) bergen2019machineBergen, K., Johnson, P., Maarten, V. Beroza, G.  2019\.  Machine learning for data-driven discovery in solid Earth geoscience Machine learning for data-driven discovery in solid earth geoscience.  Science3636433.
  * Beven Kirkby (1979) Beven1979Beven, K. Kirkby, M.  1979\.  A physically based, variable contributing area model of basin hydrology A physically based, variable contributing area model of basin hydrology.  Hydrological Sciences Journal24143–69.
  * Bookhagen (2010) Bookhagen2010Bookhagen, B.  2010\.  Appearance of extreme monsoonal rainfall events and their impact on erosion in the Himalaya Appearance of extreme monsoonal rainfall events and their impact on erosion in the Himalaya.  Geomatics, Natural Hazards and Risk1137-50.
  * Bookhagen Burbank (2010) BookhagenBurbank2010Bookhagen, B. Burbank, D.  2010\.  Toward a complete Himalayan hydrological budget: Spatiotemporal distribution of snowmelt and rainfall and their impact on river discharge Toward a complete Himalayan hydrological budget: Spatiotemporal distribution of snowmelt and rainfall and their impact on river discharge.  Journal of Geophysical Research: Earth Surface115F3.
  * Bouchard Youssef (2018) khaled_r9Bouchard, LS. Youssef, K.  201804 24.  Feature-preserving noise removal. Feature-preserving noise removal.  Google Patents.  US Patent 9,953,246
  * Burgess . (2012) Burgess2012Burgess, W., Yin, A., Dubey, C., Shen, ZK. Kelty, T.  2012\.  Holocene shortening across the Main Frontal Thrust zone in the eastern Himalaya Holocene shortening across the Main Frontal Thrust zone in the eastern Himalaya.  Earth and Planetary Science Letters357152–167.
  * Chowdhuri . (2021) chowdhuri2021Chowdhuri, I., Pal, S., Chakrabortty, R., Malik, S., Das, B. Roy, P.  2021\.  Torrential rainfall-induced landslide susceptibility assessment using machine learning and statistical methods of eastern Himalaya Torrential rainfall-induced landslide susceptibility assessment using machine learning and statistical methods of eastern himalaya.  Natural Hazards1–26.
  * Conforti . (2014) lsnn1Conforti, M., Pascale, S., Robustelli, G. Sdao, F.  2014\.  Evaluation of prediction capability of the artificial neural networks for mapping landslide susceptibility in the Turbolo River catchment (northern Calabria, Italy) Evaluation of prediction capability of the artificial neural networks for mapping landslide susceptibility in the turbolo river catchment (northern calabria, italy).  CATENA113236 - 250.
  * Coudurier-Curveur . (2020) Coudurier2020Coudurier-Curveur, A., Tapponnier, P., Okal, E., Van der Woerd, J., Kali, E., Choudhury, S.Karakaş, Ç.  2020\.  A composite rupture model for the great 1950 Assam earthquake across the cusp of the East Himalayan Syntaxis A composite rupture model for the great 1950 Assam earthquake across the cusp of the East Himalayan Syntaxis.  Earth and Planetary Science Letters531115928.
  * Cui . (2019) cui2019costCui, Y., Cheng, D., Choi, C., Jin, W., Lei, Y. Kargel, J.  2019\.  The cost of rapid and haphazard urbanization: lessons learned from the Freetown landslide disaster The cost of rapid and haphazard urbanization: lessons learned from the freetown landslide disaster.  Landslides1661167–1176.
  * Devkota . (2013) devkota2013Devkota, K., Regmi, A., Pourghasemi, H., Yoshida, K., Pradhan, B., Ryu, I.Althuwaynee, O.  2013\.  Landslide susceptibility mapping using certainty factor, index of entropy and logistic regression models in GIS and their comparison at Mugling–Narayanghat road section in Nepal Himalaya Landslide susceptibility mapping using certainty factor, index of entropy and logistic regression models in gis and their comparison at mugling–narayanghat road section in nepal himalaya.  Natural hazards651135–165.
  * Dietrich . (2001) Dietrich2001Dietrich, W., Bellugi, D. Real De Asua, R.  2001\.  Validation of the shallow landslide model, SHALSTAB, for forest management Validation of the shallow landslide model, SHALSTAB, for forest management.  Water Science and Application2195-227.
  * Dietrich . (1995) Dietrich1995Dietrich, W., Reiss, R., Hsu, M. Montgomery, D.  1995\.  A process-based model for colluvial soil depth and shallow landsliding using digital elevation data A process-based model for colluvial soil depth and shallow landsliding using digital elevation data.  Hydrological Processes93-4383-400.
  * Dikshit . (2020) dikshit2020pathwaysDikshit, A., Pradhan, B. Alamri, AM.  2020\.  Pathways and challenges of the application of artificial intelligence to geohazards modelling Pathways and challenges of the application of artificial intelligence to geohazards modelling.  Gondwana Research.
  * European Commission (2020) European2020European Commission.  2020\.  White paper on artificial intelligence–a European approach to excellence and trust White paper on artificial intelligence–a European approach to excellence and trust.
  * Fan . (2019) fan2019earthquakeFan, X., Scaringi, G., Korup, O., West, A., van Westen, C., Tanyas, H.Huang, R.  2019\.  Earthquake-induced chains of geologic hazards: Patterns, mechanisms, and impacts Earthquake-induced chains of geologic hazards: Patterns, mechanisms, and impacts.  Reviews of geophysics572421–503.
  * Fernandes Lona (2005) fernandes2005neuralFernandes, F. Lona, L.  2005\.  Neural network applications in polymerization processes Neural network applications in polymerization processes.  Brazilian Journal of Chemical Engineering223401–418.
  * Friedman (2001) friedman2001greedyFriedman, J.  2001\.  Greedy function approximation: a gradient boosting machine Greedy function approximation: a gradient boosting machine.  Annals of statistics1189–1232.
  * Froude Petley (2018) froude2018globalFroude, M. Petley, D.  2018\.  Global fatal landslide occurrence from 2004 to 2016 Global fatal landslide occurrence from 2004 to 2016.  Natural Hazards and Earth System Sciences1882161–2181.
  * Ghorbanzadeh . (2019) ghorbanzadeh2019evaluationGhorbanzadeh, O., Blaschke, T., Gholamnia, K., Meena, S., Tiede, D. Aryal, J.  2019\.  Evaluation of different machine learning methods and deep-learning convolutional neural networks for landslide detection Evaluation of different machine learning methods and deep-learning convolutional neural networks for landslide detection.  Remote Sensing112196.
  * Gómez Kavzoglu (2005) lsnn2Gómez, H. Kavzoglu, T.  2005\.  Assessment of shallow landslide susceptibility using artificial neural networks in Jabonosa River Basin, Venezuela Assessment of shallow landslide susceptibility using artificial neural networks in jabonosa river basin, venezuela.  Engineering Geology78111 - 27.  http://www.sciencedirect.com/science/article/pii/S0013795204002236
  * Gunning . (2019) refx3Gunning, D., Stefik, M., Choi, J., Miller, T., Stumpf, S. Yang, GZ.  2019\.  XAI—Explainable artificial intelligence XAI—explainable artificial intelligence.  Science Robotics437.  10.1126/scirobotics.aay7120
  * Guzzetti . (1999) guzzetti1999Guzzetti, F., Carrara, A., Cardinali, M. Reichenbach, P.  1999\.  Landslide hazard evaluation: a review of current techniques and their application in a multi-scale study, Central Italy Landslide hazard evaluation: a review of current techniques and their application in a multi-scale study, Central Italy.  Geomorphology311-4181–216.
  * Haproff . (2020) Haproff2020Haproff, P., Odlum, M., Zuza, A., Yin, A. Stockli, D.  2020\.  Structural and Thermochronologic Constraints on the Cenozoic Tectonic Development of the Northern Indo-Burma Ranges Structural and thermochronologic constraints on the Cenozoic tectonic development of the northern Indo-Burma Ranges.  Tectonics39e2020TC006231.
  * Haproff . (2019) Haproff2019Haproff, P., Zuza, A., Yin, A., Harrison, T., Manning, C., Dubey, C.Chen, J.  2019\.  Geologic framework of the northern Indo-Burma Ranges and lateral correlation of Himalayan-Tibetan lithologic units across the eastern Himalayan syntaxis Geologic framework of the northern Indo-Burma ranges and lateral correlation of Himalayan-Tibetan lithologic units across the eastern Himalayan syntaxis.  Geosphere153856-881.
  * Hastie Tibshirani (1990) hastie1990generalizedHastie, T. Tibshirani, R.  1990\.  Generalized additive models Generalized additive models ( 43).  CRC press.
  * Hinton . (2015) kdrefHinton, G., Vinyals, O. Dean, J.  2015\.  Distilling the knowledge in a neural network Distilling the knowledge in a neural network.  arXiv:1503.02531.
  * Hooker (2004) adnet1Hooker, G.  2004\.  Discovering Additive Structure in Black Box Functions Discovering additive structure in black box functions.  Proceedings of the Tenth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining Proceedings of the tenth acm sigkdd international conference on knowledge discovery and data mining ( 575–580).
  * AL. Huang Montgomery (2014) huang2014topographicHuang, AL. Montgomery, D.  2014\.  Topographic locations and size of earthquake-and typhoon-generated landslides, Tachia River, Taiwan Topographic locations and size of earthquake-and typhoon-generated landslides, tachia river, taiwan.  Earth Surface Processes and Landforms393414–418.
  * R. Huang Fan (2013) huang2013landslideHuang, R. Fan, X.  2013\.  The landslide story The landslide story.  Nature Geoscience65325–326.
  * Iverson (2000) iverson2000landslideIverson, R.  2000\.  Landslide triggering by rain infiltration Landslide triggering by rain infiltration.  Water resources research3671897–1910.
  * Kavzoglu . (2014) Kavzoglu2014Kavzoglu, T., Sahin, E. Colkesen, I.  2014\.  Landslide susceptibility mapping using GIS-based multi-criteria decision analysis, support vector machines, and logistic regression Landslide susceptibility mapping using GIS-based multi-criteria decision analysis, support vector machines, and logistic regression.  Landslides113425-439.
  * Kent Dasgupta (2004) Kent2004Kent, W. Dasgupta, U.  2004\.  Structural evolution in response to fold and thrust belt tectonics in northern Assam. A key to hydrocarbon exploration in the Jaipur anticline area Structural evolution in response to fold and thrust belt tectonics in northern Assam. a key to hydrocarbon exploration in the Jaipur anticline area.  Marine and Petroleum Geology217785-803.
  * Kingma Ba (2014) kingma2014Kingma, D. Ba, J.  2014\.  Adam: A method for stochastic optimization Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * Kirschbaum . (2020) kirschbaum2020changesKirschbaum, D., Kapnick, S., Stanley, T. Pascale, S.  2020\.  Changes in extreme precipitation and landslides over High Mountain Asia Changes in extreme precipitation and landslides over high mountain asia.  Geophysical Research Letters474e2019GL085347.
  * Krizhevsky . (2012) Krizhevsky2012Krizhevsky, A., Sutskever, I. Hinton, G.  2012\.  Imagenet classification with deep convolutional neural networks Imagenet classification with deep convolutional neural networks.  Advances in neural information processing systems Advances in neural information processing systems ( 1097-1105).
  * Larsen Montgomery (2012) Larsen2012Larsen, I. Montgomery, D.  2012\.  Landslide erosion coupled to tectonics and river incision Landslide erosion coupled to tectonics and river incision.  Nature Geoscience57468–473.
  * Larsen . (2010) Larsen2010Larsen, I., Montgomery, D. Korup, O.  2010\.  Landslide erosion controlled by hillslope material Landslide erosion controlled by hillslope material.  Nature Geoscience34247–251.
  * Lawrence . (1997) Lawrence1997Lawrence, S., Giles, C., Tsoi, A. Back, A.  1997\.  Face recognition: A convolutional neural-network approach Face recognition: A convolutional neural-network approach.  IEEE transactions on neural networks8198-113.
  * Le . (2011) khaled_r2Le, Q., Ngiam, J., Coates, A., Lahiri, A., Prochnow, B. Ng, A.  2011\.  On optimization methods for deep learning On optimization methods for deep learning.  Proceedings of the 28 th International Conference on Machine Learning, Bellevue, WA, USA Proceedings of the 28 th international conference on machine learning, bellevue, wa, usa ( 265-272).
  * Lee (2005) Lee2005Lee, S.  2005\.  Application of logistic regression model and its validation for landslide susceptibility mapping using GIS and remote sensing data Application of logistic regression model and its validation for landslide susceptibility mapping using GIS and remote sensing data.  International Journal of Remote Sensing2671477-1491.
  * Lee . (2004) lee2004determinationLee, S., Ryu, JH., Won, JS. Park, HJ.  2004\.  Determination and application of the weights for landslide susceptibility mapping using an artificial neural network Determination and application of the weights for landslide susceptibility mapping using an artificial neural network.  Engineering Geology713-4289–302.
  * Lee Sambath (2006) Lee2006Lee, S. Sambath, T.  2006\.  Landslide susceptibility mapping in the Damrei Romel area, Cambodia using frequency ratio and logistic regression models Landslide susceptibility mapping in the Damrei Romel area, Cambodia using frequency ratio and logistic regression models.  Environmental Geology506847–855.
  * Leonarduzzi . (2021) leonarduzzi2021numericalLeonarduzzi, E., Maxwell, R., Mirus, BB. Molnar, P.  2021\.  Numerical analysis of the effect of subgrid variability in a physically based hydrological model on runoff, soil moisture, and slope stability Numerical analysis of the effect of subgrid variability in a physically based hydrological model on runoff, soil moisture, and slope stability.  Water Resources Research574e2020WR027326.
  * Li . (2020) refx1Li, XH., Cao, C., Shi, Y., Bai, W., Gao, H., Qiu, L.Chen, L.  2020\.  A Survey of Data-driven and Knowledge-aware eXplainable AI A survey of data-driven and knowledge-aware eXplainable AI.  IEEE Transactions on Knowledge and Data Engineering.
  * Mandal Mandal (2018) mandal2018Mandal, S. Mandal, K.  2018\.  Modeling and mapping landslide susceptibility zones using GIS based multivariate binary logistic regression (LR) model in the Rorachu river basin of eastern Sikkim Himalaya, India Modeling and mapping landslide susceptibility zones using GIS based multivariate binary logistic regression (LR) model in the Rorachu river basin of eastern Sikkim Himalaya, India.  Modeling Earth Systems and Environment4169–88.
  * Marc Hovius (2015) marc2015amalgamationMarc, O. Hovius, N.  2015\.  Amalgamation in landslide maps: effects and automatic detection. Amalgamation in landslide maps: effects and automatic detection.  Natural Hazards & Earth System Sciences154.
  * McGuire . (2016) mcguire2016elucidatingMcGuire, L., Rengers, F., Kean, J., Coe, J., Mirus, B., Baum, R. Godt, J.  2016\.  Elucidating the role of vegetation in the initiation of rainfall-induced shallow landslides: Insights from an extreme rainfall event in the Colorado Front Range Elucidating the role of vegetation in the initiation of rainfall-induced shallow landslides: Insights from an extreme rainfall event in the Colorado front range.  Geophysical Research Letters43179084–9092.
  * Meunier . (20081) Menuier2008Meunier, P., Hovius, N. Haines, J.  20081\.  Topographic site effects and the location of earthquake induced landslides Topographic site effects and the location of earthquake induced landslides.  Earth and Planetary Science Letters2753-4221-232.
  * Meunier . (20082) meunier2008topographicMeunier, P., Hovius, N. Haines, J.  20082\.  Topographic site effects and the location of earthquake induced landslides Topographic site effects and the location of earthquake induced landslides.  Earth and Planetary Science Letters2753-4221–232.
  * Montavon . (2012) khaled_r5Montavon, G., Orr, G. Müller, KR.  2012\.  Neural Networks: Tricks of the Trade Neural Networks: Tricks of the Trade ( 7700).  Springer.
  * Montgomery Dietrich (1994) Montgomery1994Montgomery, D. Dietrich, W.  1994\.  A physically based model for the topographic control on shallow landsliding A physically based model for the topographic control on shallow landsliding.  Water Resources Research3041153-1171.
  * Montgomery . (1998) Montgomery1998Montgomery, D., Sullivan, K. Greenberg, H.  1998\.  Regional test of a model for shallow landsliding Regional test of a model for shallow landsliding.  Hydrological processes126943–955.
  * Moon . (2011) Moon2011Moon, S., Chamberlain, C., Blisniuk, K., Levine, N., Rood, D. Hilley, G.  2011\.  Climatic control of denudation in the deglaciated landscape of the Washington Cascades Climatic control of denudation in the deglaciated landscape of the Washington Cascades.  Nature Geoscience47469-473.
  * Orland . (2020) Orland2020Orland, E., Roering, J., Thomas, M. Mirus, B.  2020\.  Deep Learning as a Tool to Forecast Hydrologic Response for Landslide-Prone Hillslopes Deep learning as a tool to forecast hydrologic response for landslide-prone hillslopes.  Geophysical Research Letters4716e2020GL088731.  10.1029/2020GL088731
  * Ozyildirim Kiran (2020) khaled_r1Ozyildirim, B. Kiran, M.  2020\.  Do optimization methods in deep learning applications matter? Do optimization methods in deep learning applications matter?  arXiv:2002.12642.
  * Parker . (2011) parker2011massParker, R., Densmore, A., Rosser, N., De Michele, M., Li, Y., Huang, R.Petley, D.  2011\.  Mass wasting triggered by the 2008 Wenchuan earthquake is greater than orogenic growth Mass wasting triggered by the 2008 wenchuan earthquake is greater than orogenic growth.  Nature Geoscience47449–452.
  * Petley (2012) petley2012globalPetley, D.  2012\.  Global patterns of loss of life from landslides Global patterns of loss of life from landslides.  Geology4010927–930.
  * Phong . (2019) phong2019landslidePhong, T., Phan, T., Prakash, I., Singh, S., Shirzadi, A., Chapi, K.Pham, B.  2019\.  Landslide susceptibility modeling using different artificial intelligence methods: A case study at Muong Lay district, Vietnam Landslide susceptibility modeling using different artificial intelligence methods: A case study at Muong Lay district, Vietnam.  Geocarto International1–24.
  * Pradhan (2013) Pradhan2013Pradhan, B.  2013\.  A comparative study on the predictive ability of the decision tree, support vector machine and neuro-fuzzy models in landslide susceptibility mapping using GIS A comparative study on the predictive ability of the decision tree, support vector machine and neuro-fuzzy models in landslide susceptibility mapping using GIS.  Computers & Geosciences51350-365.
  * Prakash . (2020) prakash2020mappingPrakash, N., Manconi, A. Loew, S.  2020\.  Mapping landslides on EO data: Performance of deep learning models vs. traditional machine learning models Mapping landslides on eo data: Performance of deep learning models vs. traditional machine learning models.  Remote Sensing123346.
  * Radbruch-Hall (1982) radbruch1982landslideRadbruch-Hall, DH.  1982\.  Landslide overview map of the conterminous United States Landslide overview map of the conterminous united states ( 1183).  US Government Printing Office.
  * Regmi . (2014) Regmi2014Regmi, A., Devkota, K., Yoshida, K., Pradhan, B., Pourghasemi, H., Kumamoto, T. Akgun, A.  2014\.  Application of frequency ratio, statistical index, and weights-of-evidence models and their comparison in landslide susceptibility mapping in Central Nepal Himalaya Application of frequency ratio, statistical index, and weights-of-evidence models and their comparison in landslide susceptibility mapping in central Nepal Himalaya.  Arabian Journal of Geosciences72725-742.
  * Reichenbach . (2018) reichenbach2018reviewReichenbach, P., Rossi, M., Malamud, B., Mihir, M. Guzzetti, F.  2018\.  A review of statistically-based landslide susceptibility models A review of statistically-based landslide susceptibility models.  Earth-Science Reviews18060–91.
  * Rudin (2019) rudin2019stopRudin, C.  2019\.  Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead.  Nature Machine Intelligence15206–215.
  * Salvi . (2020) Salvi2020Salvi, D., Mathew, G., Kohn, B., Pande, K. Borgohain, B.  2020\.  Thermochronological insights into the thermotectonic evolution of Mishmi Hills across the Dibang Valley, NE Himalayan Syntaxis Thermochronological insights into the thermotectonic evolution of Mishmi Hills across the Dibang Valley, NE Himalayan Syntaxis.  Journal of Asian Earth Sciences190104158.
  * Scherler . (2014) Scherler2014Scherler, D., Bookhagen, B. Strecker, M.  2014\.  Tectonic control on 10Be‐derived erosion rates in the Garhwal Himalaya, India Tectonic control on 10Be‐derived erosion rates in the Garhwal Himalaya, India.  Journal of Geophysical Research: Earth Surface119283-105.
  * Schmidt . (2003) Schmidt2003Schmidt, J., Evans, I. Brinkmann, J.  2003\.  Comparison of polynomial models for land surface curvature calculation Comparison of polynomial models for land surface curvature calculation.  International Journal of Geographical Information Science178797-814.
  * Schwanghart Scherler (2014) Schwanghart2014Schwanghart, W. Scherler, D.  2014\.  TopoToolbox 2–MATLAB-based software for topographic analysis and modeling in Earth surface sciences TopoToolbox 2–MATLAB-based software for topographic analysis and modeling in Earth surface sciences.  Earth Surface Dynamics211-7.
  * Stanley Kirschbaum (2017) stanley2017heuristicStanley, T. Kirschbaum, DB.  2017\.  A heuristic approach to global landslide susceptibility mapping A heuristic approach to global landslide susceptibility mapping.  Natural hazards871145–164.
  * Stine (1995) Stine1995Stine, R.  1995\.  Graphical interpretation of variance inflation factors Graphical interpretation of variance inflation factors.  The American Statistician49153-56.
  * Stock Dietrich (2003) stock2003valleyStock, J. Dietrich, W.  2003\.  Valley incision by debris flows: Evidence of a topographic signature Valley incision by debris flows: Evidence of a topographic signature.  Water Resources Research394.
  * Stock Dietrich (2006) stock2006erosionStock, J. Dietrich, W.  2006\.  Erosion of steepland valleys by debris flows Erosion of steepland valleys by debris flows.  Geological Society of America Bulletin1189-101125–1148.
  * H. Tan Lim (2019) khaled_r4Tan, H. Lim, K.  2019jun.  Review of second-order optimization techniques in artificial neural networks backpropagation Review of second-order optimization techniques in artificial neural networks backpropagation.  IOP Conference Series: Materials Science and Engineering495012003.
  * S. Tan . (2018) adnet2Tan, S., Caruana, R., Hooker, G., Koch, P. Gordo, A.  2018\.  Learning global additive explanations for neural nets using model distillation Learning global additive explanations for neural nets using model distillation.  arXiv:1801.08640.
  * Taylor Yin (2009) Taylor2009Taylor, M. Yin, A.  2009\.  Active structures of the Himalayan-Tibetan orogen and their relationships to earthquake distribution, contemporary strain field, and Cenozoic volcanism Active structures on the Tibetan Plateau and surrounding regions Active structures of the Himalayan-Tibetan orogen and their relationships to earthquake distribution, contemporary strain field, and Cenozoic volcanism active structures on the Tibetan Plateau and surrounding regions.  Geosphere53199-214.
  * Tien Bui . (2012) bui2012landslideTien Bui, D., Pradhan, B., Lofman, O., Revhaug, I. Dick, O.  2012\.  Landslide susceptibility assessment in the Hoa Binh province of Vietnam: a comparison of the Levenberg–Marquardt and Bayesian regularized neural networks Landslide susceptibility assessment in the hoa binh province of vietnam: a comparison of the levenberg–marquardt and bayesian regularized neural networks.  Geomorphology17112–29.
  * Tien Bui . (2019) tien2019shallowTien Bui, D., Shahabi, H., Omidvar, E., Shirzadi, A., Geertsema, M., Clague, J.Lee, S.  2019\.  Shallow landslide prediction using a novel hybrid functional machine learning algorithm Shallow landslide prediction using a novel hybrid functional machine learning algorithm.  Remote Sensing118931.
  * Tien Bui . (2020) bui2020comparingTien Bui, D., Tsangaratos, P., Nguyen, VT., Van Liem, N. Trinh, P.  2020\.  Comparing the prediction performance of a Deep Learning Neural Network model with conventional machine learning models in landslide susceptibility assessment Comparing the prediction performance of a deep learning neural network model with conventional machine learning models in landslide susceptibility assessment.  Catena188104426.
  * Toms . (2020) toms2020physicallyToms, B., Barnes, E. Ebert-Uphoff, I.  2020\.  Physically interpretable neural networks for the geosciences: Applications to earth system variability Physically interpretable neural networks for the geosciences: Applications to earth system variability.  Journal of Advances in Modeling Earth Systems129e2019MS002002.
  * Tucker . (1986) Tucker1986Tucker, C., Fung, I., Keeling, C. Gammon, R.  1986\.  Relationship between atmospheric CO2 variations and a satellite-derived vegetation index Relationship between atmospheric CO2 variations and a satellite-derived vegetation index.  Nature3196050195–199.
  * United States Geological Survey EarthExplorer (2020) EarthExplorerUnited States Geological Survey EarthExplorer.  2020\.  accessed September 26, 2020, https://earthexplorer.usgs.gov/. accessed September 26, 2020, https://earthexplorer.usgs.gov/.