Uncertainty quantification of landslide runout motion considering soil interdependent anisotropy and fabric orientation

Abstract

Natural soils often exhibit an anisotropic fabric pattern as a result of soil deposition, weathering, or filling. This paper aims to investigate the effect of soil interdependent anisotropy and fabric orientation on runout motions of landslides and evaluate the most critical fabric orientation for the post-failure behavior. The shear strength properties of soil deposit (i.e., cohesion $c$ and friction angle $\varphi$) are modeled as negatively cross correlated bivariate random fields. The results reveal that the spatial variability and the negative cross-correlation of $c$ and $\varphi$ notably influence the post-failure behavior. In addition, the rotation of soil layer orientation significantly affects the runout motion. Based on the analyses, the deposition orientation of ${30}^{^\circ }$ is identified to produce the highest mean value and standard deviation of the runout distance. The findings from this study highlight the importance of considering the orientation of soil stratification, rather than only the magnitude of shear strength, in assessing the post-failure behavior of a landslide.

Introduction

A large number of geohazards, such as landslides, occur every year in mountainous areas in China, Japan, Italy, and other places in the world. In these hazards, hydro-mechanical characteristics of heterogeneous natural soils exhibit considerable spatial variability (Nezhad et al. 2011), even within homogeneous soil layers, which can have considerable impacts on instability of slopes as well as their failure likelihood, mechanisms, and post-failure runout motions (Zhu and Zhang 2013; Zhu et al. 2019). Diverse failure consequences can be observed in a varying soil slope consisting of weak materials (e.g., weak intercalation layer) because the failure surface tends to seek out the weakest path. In many occasions, due to extensive runout motions of landslides during the post-failure stage, substantial destruction can be observed on nearby structures (Yin et al. 2009, 2016), which leads to severe threat and significant damages to human life and properties (Ma et al. 2018; He et al. 2019; Wang et al. 2021). Studying the post-failure responses in heterogeneous soil slopes is rather complicated due to large deformations and interdependencies in spatially varying shear strengths. Therefore, it is necessary to propose a reliable predictive model to assess the probable post-failure consequences, especially the landslide runout distance which is of great concern and employ it for timely planning of disaster mitigation measures.

In recent years, probabilistic studies primarily aim to describe the slope instability or corresponding probability of failure (${P}_{f}$) using limit equilibrium method (LEM) (e.g., Griffiths and Fenton 2004; Jiang et al. 2014) or finite element method (FEM) integrated with random field (RF) modeling (e.g., Nezhad et al. 2018a, 2018b, 2019). It has been recognized that consideration of natural heterogeneity in soil properties has a significant effect on the analysis of soil deformation. However, since abovementioned modeling methods, for assessment of slope instability, generally stop at the pre-failure stage, limited research can be found on probabilistic analyses of landslides’ post-failure behavior. The previous studies only provide limited information about the consequences of landslides in spatially varying soils and do not quantify their possible secondary motions. The main reason is related to the numerical instability of the classical Galerkin FEM modeling algorithms due to mesh distortion or twisting that occur when simulating large deformation problems (Hardcastle et al. 2019). Recent progress in developing particle-based methods, such as material point method (MPM) (e.g., Yerro, 2015), made it possible to tackle the numerical instability issues associated with the mesh distortion problems. The MPM has been widely adopted to simulate landslides, slope failures, embankment collapses, and other large strain geotechnical problems (e.g., Andersen and Andersen 2010; Llano-Serna et al. 2016; Li et al. 2016). However, in the previous MPM studies, the materials have been mainly assumed to be simply homogeneous and isotropic, and the geomaterials’ responses have been modeled deterministically using constant values of soil parameters. Recently, in few studies random MPM was used to investigate the slope failure process, for example Liu et al. (2019) utilized the LEM and MPM combined with Monte-Carlo simulations (MCS) to quantify the failure probability of landslides and predicted the slope failure modes. However, these models only considered one shear strength parameter (e.g., undrained shear strength) as a random variable, and as such they are unable to deal with deposits with frictional characteristics.

In reality, interdependent cohesion and friction angle of soils vary significantly in space, which are the significant factors influencing the deformation behavior (Wang et al. 2020a). Especially, the cross-correlation between the shear strength parameters on geo-structures cannot be ignored (Wu 2015; Wang et al. 2020b), which has significant influence on soil deformation. Cohesion, friction angle, and their cross-correlation effect on post-failure behavior of landslides have not been previously investigated in detail. Furthermore, previous probabilistic studies of slope failures or landslides mainly focused on soils with isotropic or horizontally deposited strata (Liu et al. 2019); whereas, in nature slopes with tilted stratifications are often observed (Ma et al. 2018; He et al. 2021). The geometric relationship between a slope and underlying layers/strata can considerably influence the landslide behavior. Generally, from the perspective of soil profile, there are two main types of slopes’ bedding that can categorize their topographical features: (i) a dip slope (also termed as cataclinal slope) as shown in Fig. 1a, (ii) an anti-dip slope (also termed as anaclinal slope or reverse slope) as shown in Fig. 1b. In engineering geology, it is commonly admitted that the dip slope is associated with planar sliding failure, which is favorable geological structure to form large-scale consequent bedding slides, while the anti-dip slope generally fails in the form of toppling. It should be noted that the failure probability and kinematic likelihood of failures in these slope types also vary. It has been found that the reverses slopes with the dip direction of the strata being against the dip direction of the slope usually have better stability or reliability compared to when the dip direction of the strata is the same as the dip direction of the slope. Recently, the effect of anisotropy and fabric orientation of soil layers on the probabilistic stability or reliability analysis of slope failure have been considered in a number of works. For example, Griffiths et al. (2009) found that the probability of failure is high when the dipping direction of soil layers is aligned parallel to the slope surface. Zhu and Zhang (2013) discussed the typical patterns of spatial variability of geomaterials, such as isotropy, transverse anisotropy, rotated anisotropy, and general anisotropy. Based on the analysis of a drained slope with fabric orientation of different angles, they found that the rotated anisotropy in soil properties can significantly influence the probability of failure as well as the failure mechanisms (Zhu et al. 2019). Cheng et al. (2018) adopted the random finite difference method to estimate failure consequence of slopes with different rotational angles. Despite these previous works, the influence mechanisms of soil anisotropy and fabric orientation on landslides and their post-failure behavior are still unclear.

Fig. 1

Two main types of bedding slopes. a dip slope, b anti-dip slope

This paper proposes the modeling framework for the probabilistic analysis of post-failure behavior of landslides with interdependent anisotropy and fabric orientation in cohesive-frictional soil, where the cross-correlated bivariate RFs (interdependent cohesion and friction angle parameters) are incorporated into the generalized interpolation MPM (GIMP) through an MCS algorithm. The layout of the paper is as follows. The “Methodology” section describes the probabilistic framework in which the cross-correlated bivariate RFs are generated for different auto-correlation structures via Cholesky matrix decomposition method coupled with Latin hypercube sampling method. Then, the “Illustrative example” section describes a cohesive-frictional soil slope that is used as the main example to demonstrate the performance and validity of the proposed framework. Thereafter, “Deterministic analysis of landslide runout motion” is discussed for the example landslide, and “Comparative analyses and stochastic assessment” are carried out for post-failure motions under different parametric settings. In the section on “Effect of auto-correlation structures,” sensitivity analyses are performed on three auto-correlation structures of soil shear strength parameters to investigate how an anisotropically deposited strata influences the runout motions, and to identify the critical auto-correlation structure which leads to the most extensive runout motion behavior of heterogeneous landslides. The corresponding statistic characteristics of the runout distances, including their means, variance, and corresponding best-fit marginal distributions, are eventually obtained, and discussed.

Methodology

Generation of bivariate random fields

Auto-correlation structures

RF modeling has been commonly adopted to mathematically characterize the spatial variability of shear strength in soils. In the RF, an important measure of the variability is the auto-correlation distance, which may not be constant for different directions in a soil mass. The spatial variability of soils’ shear strength could entail auto-correlation distance ranging from less than 1 m to tens of meters along different directions, which produces diverse auto-correlation structures. Figure 2 shows the six typical patterns of spatially varying geomaterial profiles (Zhu and Zhang, 2013; Cheng et al. 2018; Zhu et al. 2019). Figure 2a is an isotropic case, where the material properties are independent of the direction. In Fig. 2b, the distribution of soil’s structure is directional along two orthogonal principal directions, and the case is defined as transverse anisotropy. The longer line with an arrow represents the direction with less variations of soil properties (i.e., horizontal homogeneity), and the shorter one represents the rapid variations of structure (i.e., vertical heterogeneity). This type of soil structure is commonly found in geological horizontal depositions, fluvial processes, and soil compaction. Figure 2c presents the rotated anisotropy case caused by geological processes, in which the two principal directions are rotated by an angle $\beta$ to the horizon. In Fig. 2d, the two principal directions form an angle $\eta$ which is not orthogonal; this case is defined as a general anisotropy case. Figure 2e and f show the rotated general anisotropy and a combination case involving two of the above structures in a profile, respectively. As general anisotropy and general rotated anisotropy cases are less likely to exist in a soil mass, this work only focuses on the influence of the first three auto-correlation structures in landslides.

Fig. 2

Illustration of different auto-correlation structures. a isotropy, b transverse anisotropy, c rotated anisotropy, d general anisotropy, e general rotated anisotropy, f combination

In geotechnical RFs, two theoretical expressions are commonly employed to describe the auto-correlation structure of soil mass in two dimensions: single exponential ACFs (Li et al. 2015) and squared exponential ACFs (Cheng et al. 2018; Masoudian et al. 2019). According to previous studies (e.g., Cheng et al. 2018), the single exponential ACFs do not satisfy the requirement of transverse anisotropy and rotated anisotropy structures; therefore, in this work the squared exponential ACFs are used. The expressions of the squared exponential ACFs correlation coefficient $\rho ({\tau }_{x},{\tau }_{y})$ and auto-correlation distance ${\theta }_{\phi }$ are given as.

Isotropy:

$$\begin{array}{c}\rho \left({\tau }_{x},{\tau }_{y}\right)=exp\left[-\frac{\sqrt{{\tau }_{x}^{2}+{\tau }_{y}^{2}}}{\theta }\right]\end{array}$$

(1)

$$\begin{array}{c}{\theta }_{\phi }={\theta }_{1}={\theta }_{2}\end{array}$$

(2)

Transverse anisotropy:

$$\begin{array}{c}\rho \left({\tau }_{x},{\tau }_{y}\right)=exp\left[-\sqrt{{\left(\frac{{\tau }_{x}}{{\theta }_{1}}\right)}^{2}+{\left(\frac{{\tau }_{y}}{{\theta }_{2}}\right)}^{2}}\right]\end{array}$$

(3)

$$\begin{array}{c}{\theta }_{\phi }=\sqrt{\frac{{\theta }_{1}^{2}{\theta }_{2}^{2}\left(1+{tan}^{2}\phi \right)}{{\theta }_{2}^{2}+{\theta }_{1}^{2}{tan}^{2}\phi }}\end{array}$$

(4)

Rotated anisotropy:

$$\begin{array}{c}\rho \left({\tau }_{x},{\tau }_{y}\right)=exp\left[-\sqrt{{\left(\frac{{\tau }_{x}cos\beta +{\tau }_{y}cos\beta }{{\theta }_{1}}\right)}^{2}+{\left(\frac{-{\tau }_{x}cos\beta +{\tau }_{y}cos\beta }{{\theta }_{2}}\right)}^{2}}\right]\end{array}$$

(5)

$$\begin{array}{c}{\theta }_{\phi }=\sqrt{\frac{{\theta }_{1}^{2}{\theta }_{2}^{2}\left(1+{tan}^{2}\phi \right)}{{\theta }_{2}^{2}{(cos\beta +sin\beta tan\phi )}^{2}+{{\theta }_{1}^{2}(-sin\beta +cos\beta tan\phi )}^{2}}}\end{array}$$

(6)

where $\rho (\tau )$ is the auto-correlation coefficients, and ${\tau }_{x}$ and ${\tau }_{y}$ are the relative distances between any two points in horizontal and vertical directions, respectively; ${\theta }_{1}$ and ${\theta }_{2}$ are the auto-correlation distance in principal directions;$\phi$ is the directional angle; $\beta$ is the rotated angle of the principal axis in the rotated anisotropy structure, when $\beta ={0}^{^\circ } \mathrm{or} {180}^{^\circ }$ denotes the tranverse anisotropy structure.

Once the ACF is determined, the domain $\Omega$ can be discretized into $n$ elements of the RF with the centroid coordinates of the elements specified by $({x}_{i},{y}_{i})$. The auto-correlation matrix ${{\varvec{C}}}_{n\times n}$, representing the spatial variability of the properties (Nezhad, 2010), can be obtained with respect to the elements as

$${{\varvec{C}}}_{n\times n}=\left[\begin{array}{cccc}1& \rho ({\tau }_{{x}_{12}},{\tau }_{{y}_{12}})& \dots & \rho ({\tau }_{{x}_{1n}},{\tau }_{{y}_{1n}})\\ \rho ({\tau }_{{x}_{21}},{\tau }_{{y}_{21}})& 1& \dots & \rho ({\tau }_{{x}_{2n}},{\tau }_{{y}_{2n}})\\ \vdots & \vdots & \ddots & \vdots \\ \rho ({\tau }_{{x}_{n1}},{\tau }_{{y}_{n1}})& \rho ({\tau }_{{x}_{n2}},{\tau }_{{y}_{n2}})& \dots & 1\end{array}\right]$$

(7)

Subsequently, the ${{\varvec{C}}}_{n\times n}$ can be decomposed by Cholesky matrix decomposition method (Jiang et al. 2014; Liu et al. 2017) into the product of a lower triangular matrix ${\varvec{L}}$ and a conjugate transpose ${{\varvec{L}}}^{\mathbf{T}}$ as

$$\begin{array}{c}{\mathbf{L}}\bullet {\mathbf{L}}^{\mathbf{T}}={{\varvec{C}}}_{n\times n}\end{array}$$

(8)

where $\mathbf{L}$ is the lower triangular matrix with dimension of $n\times n$. To improve the robustness and accuracy of MCS, a stratified sampling scheme, Latin hypercube sampling (LHS) method is applied in this work, which ensures that the region with prescribed standard normal distribution between 0 and 1 is uniformly divided into n non-overlapping intervals for each random variable

$$\begin{array}{c}{\xi }_{i}=\frac{\xi }{n}+\frac{i-1}{n}\end{array}$$

(9)

where $i=1, 2, 3,..., n$, and $\xi$ is a random number in the range of [0, 1]; ${\xi }_{i}$ is the random number in the ith interval. In addition, the ${\xi }_{i}$ meets the following relationship

$$\begin{array}{c}\frac{i-1}{n}<{\xi }_{i}<\frac{1}{n}\end{array}$$

(10)

where $\frac{i-1}{n}$ and $\frac{1}{n}$ are the lower bound and upper bound for the ith interval. Subsequently, a standard bivariate independent Gaussian random field ${{\varvec{X}}}^{IG}=[{\mathbf{X}}_{c}^{IG}, {\mathbf{X}}_{\varphi }^{IG}]$ of $c$ and $\varphi$ is derived as

$$\begin{array}{c}{{\varvec{X}}}^{IG}\left(x,y\right)=\mathbf{L}\bullet {\xi }_{i} \left(i=\mathrm{1,2},\dots ,N\right)\end{array}$$

(11)

where $\xi_i$ is a $n\times 2$ LHS sample matrix, and $i$ is the number of generated realizations of bivariate Gaussian RFs of $c$ and $\boldsymbol{\varphi }$.

Cross-correlation of bivariate random fields

In general, many soil properties should be involved in landslide failure analysis, which means that more than one RF is generated in a soil layer. Cohesion, $c$, and friction angle, $\varphi$, of soils are the significant factors influencing slope stability, and they are found to be negatively cross-correlated with each other with a Pearson cross-correlation coefficient ${\rho }_{c, \varphi }$. The cross-correlation matrix ${{\varvec{R}}}_{2\times 2}$ between $c$ and $\varphi$ is formed as

$$\begin{array}{c}{{\varvec{R}}}_{2\times 2}=\left[\begin{array}{cc}1& {\rho }_{c,\varphi }\\ {\rho }_{c,\varphi }& 1\end{array}\right]\end{array}$$

(12)

The cross-correlation matrix ${{\varvec{R}}}_{2\times 2}$ is decomposed by using Cholesky decomposition techniques as

$$\begin{array}{c}{{{\varvec{R}}}_{2\times 2}=\mathbf{L}}_{\rho }\cdot {\mathbf{L}}_{\rho }^{\mathbf{T}}\end{array}$$

(13)

where ${\mathbf{L}}_{\rho }$ is the lower triangular matrix of the correlation matrix ${{\varvec{R}}}_{2\times 2}$. The cross-correlated Gaussian RFs ${\mathbf{X}}^{CG}=[{\mathbf{X}}_{c}^{CG},{\mathbf{X}}_{\varphi }^{CG}]$ can be written as

$$\begin{array}{c}{\mathbf{X}}^{CG}=\left[\begin{array}{c}{\mathbf{X}}_{c}^{CG}\\ {\mathbf{X}}_{\varphi }^{CG}\end{array}\right]={\mathbf{L}}_{\rho }\cdot \left[\begin{array}{c}{\mathbf{X}}_{c}^{IG}\\ {\mathbf{X}}_{\varphi }^{IG}\end{array}\right]=\left[\begin{array}{cc}1& 0\\ {\rho }_{c,\varphi }& \sqrt{1-{\rho }_{c,\varphi }^{2}}\end{array}\right]\cdot \left[\begin{array}{c}{\mathbf{X}}_{c}^{IG}\\ {\mathbf{X}}_{\varphi }^{IG}\end{array}\right]\end{array}$$

(14)

Then, cross-correlated standard uniform RFs ${\mathbf{X}}^{CU}$ can be expressed as

$$\begin{array}{c}{\mathbf{X}}^{CU}=\Phi \left({\mathbf{X}}^{CG}\right)\end{array}$$

(15)

where $\Phi (\bullet )$ is the standard Gaussian cumulative density function (CDF). The $\Phi (\bullet )$ can rearrange the standard Gaussian RF $\mathbf{X}$ into the uniform distribution $\mathbf{U}$ within $(\mathrm{0,1})$. The cross-correlated standard uniform RFs of $c$ and $\varphi$, ${\mathbf{X}}^{CU}=[{\mathbf{X}}_{c}^{CU},{\mathbf{X}}_{\varphi }^{CU}]$, are obtained. Subsequently, the non-Gaussian cross-correlated bivariate RF ${\mathbf{X}}^{NG}$ of desired values are obtained by isoprobabilistic transformation method as

$$\begin{array}{c}{\mathbf{X}}^{NG}\left(x,y\right)={F}^{-1}\left[{\mathbf{X}}^{CU} \left(x,y\right)\right]\end{array}$$

(16)

where the ${F}^{-1}(\bullet )$ is the inverse function of its corresponding marginal cumulative distribution of the desired shear strength. If the shear strength parameters ($c$ and $\varphi$) are considered to be lognormally distributed, the realization of approximate cross-correlated lognormal RF can be obtained by exponential function

$$\begin{array}{c}{\mathbf{X}}^{LN}=\left[\begin{array}{c}{\mathbf{X}}_{c}^{LN}\\ {\mathbf{X}}_{\varphi }^{LN}\end{array}\right]=\left[\begin{array}{c}exp\left({\mu }_{ln c}+{\sigma }_{ln c}\cdot {\mathbf{X}}_{c}^{CG}\right)\\ exp\left({\mu }_{ln \varphi }+{\sigma }_{ln \varphi }\cdot {\mathbf{X}}_{\varphi }^{CG}\right)\end{array}\right]\end{array}$$

(17)

where ${\mu }_{ln c}$ and ${\sigma }_{ln c}$ are the mean and standard deviation of Gaussian random variable $ln c$, respectively; ${\mu }_{ln \varphi }$ and ${\sigma }_{ln \varphi }$ are the mean and standard deviation of Gaussian random variable $ln \varphi$, respectively. The procedure can be repeated $N$ times to obtain $N$ simulations of the bivariate RFs.

Calculation of the runout distance by GIMP

After obtaining RFs of the soil parameters, the runout distance of each bivariate RF sample can be calculated by MPM analysis. However, the original MPM has grid-crossing instability, which is caused by the discontinuous gradient of shape functions (Bardenhagen, 2002). A sudden change of the stress can be found when a material point crosses to a new cell. This deficiency can be alleviated by using generalized interpolation material point method (GIMP) that introduces alternative grid shape function, ${S}_{I}$, and particle characteristic function, ${\chi }_{p}(x)$ (Bardenhagen and Kober 2004). In this paper, the open-source GIMP program named MPM3D, from the Computational Dynamics Laboratory of Tsinghua University (Zhang et al. 2016), is used for the analysis. The implementation of computational cycle of GIMP can be found in many studies, such as Li et al. (2016, 2020), and is summarized in the Appendix section here for the sake of completeness. Meanwhile, as for the time when the stress is updated, the modified update stress last (MUSL) is adopted in this work, given its better computational stability (Nairn 2003) compared with update stress first (Bardenhagen 2002) and update stress last (Sulsky et al. 1994).

Applied soil constitutive model

The soil properties were modeled as elastic-perfectly plastic materials with the Mohr–Coulomb failure criterion. It is admitted that the strength properties are significantly degraded in mobilized soil mass during the sliding. In this work, the strain-softening behavior induced by increasing deviatoric plastic strain is used, and the corresponding softening rules are defined as

$$\begin{array}{c}c={c}_{r}+\left({c}_{p}-{c}_{r}\right){e}^{-\eta <{\varepsilon }_{ep}-{\varepsilon }_{pp}>}\end{array}$$

(18)

$$\begin{array}{c}\varphi ={\varphi }_{r}+\left({\varphi }_{p}-{\varphi }_{r}\right){e}^{-\eta <{\varepsilon }_{ep}-{\varepsilon }_{pp}>}\end{array}$$

(19)

where ${c}_{r}$ and ${c}_{p}$ are residual cohesion and peak cohesion, ${\varphi }_{r}$ and ${\varphi }_{p}$ are residual friction angle and peak friction angle, and $\eta$ and ${\varepsilon }_{ep}$ are shape factor and deviatoric plastic strain, respectively. Moreover, ${\varepsilon }_{pp}$ shows the threshold strain of softening behavior. Regarding Eqs. (18) and (19), the cohesion and friction angle of the soil remain constant for ${\varepsilon }_{ep}<{\varepsilon }_{pp}$. Once the plastic deviatoric strain exceeds the threshold strain (${\varepsilon }_{ep}>{\varepsilon }_{pp}$), the strength parameters start to reduce and tend to their residual values. Detailed description of the constitutive model in MPM can be found in Yerro et al. (2016).

Implementation and workflow

To model the probabilistic post-failure behavior of landslides that accounts for soil anisotropy, in this work, GIMP is coupled with the bivariate RF theory under MCS framework. Firstly, a RF mesh is generated for a sampling process. The CMD method with LHS is employed and a computational algorithm is written using PYTHON 3.9 to obtain bivariate cross-correlated non-Gaussian RF samples of the shear strength parameters over the RF mesh. Then, the samples are mapped onto the material points in the GIMP, in which each particle has a unique value of soil properties ($c$ and $\varphi$). The transmission process is based on the spatial relationship between the material point and the RF mesh element (e.g., a position-to-position mapping process), which is similar to the methodology used in the random FEM (Gironacci et al. 2018; Nezhad et al. 2018a, 2018b; Huang et al. 2020). Therefore, the sample values of each element in the RF mesh are assigned to the corresponding material points according to their spatial coordinates. It should be noted that each material point carries random values based on the RF, namely, $c$ and $\varphi$, repectively.

Figure 3 shows the flowchart of the proposed framework. The details are summarized in the following steps:

• 1. Identify the input information for probabilistic analysis of post-failure behavior, including slope geometry and statistics of $c$ and $\varphi$ (i.e., mean, COV, distribution type, auto-correlation structure, auto-correlation distance) and establish the GIMP deterministic slope model.

• 2. Set the total number of bivariate RFs, N; the current simulation starts from i = 1.

• 3. Discretize the domain into n elements and extract the centroid coordinates of each element. According to these coordinates, the auto-correlation matrix ${{\varvec{C}}}_{n\times n}$ can be obtained. Based on prescribed statistical information of the soil properties, one isotropic bivariate Gaussian RF with LHS ${{\varvec{X}}}_{i}^{IG}=[{{\varvec{X}}}_{c}^{IG}, {{\varvec{X}}}_{\varphi }^{IG}]$ is generated by Eq. (11) and subsequently it is transformed into a lognormal RF ${{\varvec{X}}}_{i}^{LN}\left(x,y\right)$ of $c$ and $\varphi$ by Eq. (17).

• 4. Determine the runout distance by performing a deterministic GIMP analysis of a landslide with the ith RF generated in step 3.

• 5. Update the simulation index, i, and repeat N times. It should be noted that each simulation index i would be updated and checked; if i $>$ N, abort the MC simulation; otherwise continue.

• 6. All calculated results (e.g., runout distances) are collected and post-processed to estimate the corresponding best-fit marginal distribution and statistical characteristics. The stochastic data are compared with the deterministic analysis to investigate the effect of spatial variability.

Fig. 3

Flowchart of the probabilistic analysis of landslides considering bivariate RFs

Illustrative example

An example cohesive-frictional soil slope is considered to show the feasibility of the proposed probabilistic post-failure analysis of landslides. Figure 4 shows the geometry of the slope with 15 m height and a slope gradient of 45 ${}^{\circ }$. A mesh with four-noded square elements of size 0.2 m is used, where in each element there are 4 equally spaced material points with a size of 0.1 m. The total number of material points is 10,090. The RF is generated with a cell size equal to material point domain (0.1 m), where the properties of a single cell correspond to one material point. Therefore, as RF cell is 4 times smaller than the smallest auto-correlation distance, the spatial variability is adequately captured (Nezhad et al. 2018a). The bottom of the slope is a fixed boundary, both lateral sides are roller type symmetrical boundaries, and the top boundary is free.

Fig. 4

Geometry and material point model for the example slope

As the shear strength of soil deposit is spatially distributed, the cohesion and friction angle parameters are considered cross-correlated bivariate RFs. Both point statistics of $c$ and $\varphi$ are assumed to be lognormally distributed variables to avoid negative quantities in the sample values. The COV values for $c$ and $\varphi$ are changed from 0.1 to 0.5 according to the typical values summarized by Phoon and Kulhawy (1996) and Wu (2013). For incorporating the dependence relationship between $c$ and $\varphi$, a cross-correlation coefficient ${\rho }_{c, \varphi }$ is needed. Based on the literature (Cho 2010; Wang and Akeju 2016; Fang et al. 2020), a negative cross-correlation coefficient that ranges between − 0.7 and − 0.2 is selected in this study. Meanwhile, an additional case with ${\rho }_{c, \varphi }=0$ is also considered to reveal the effect on runout distance by uncorrelated variables. All statistical moments for the soil parameters are shown in Table 1. Typical values are assigned to other parameters, i.e., Young’s modulus $E= 100 \mathrm{MPa}$, Poisson’s ratio $\nu =0.35$, and unit weight $\gamma = 20 \mathrm{kN}/{\mathrm{m}}^{3}$, as their contributions on landslide runout distance are not significant (Cheuk et al. 2013). The residual cohesion is set as 10 kPa and residual friction angle is set equal to $5{}^{\circ }$. The initial in situ stresses are generated using gravitational loading. The time increment is set equal to 7.5 $\times$ 10⁻⁴ s, and the total duration for the calculation is 10 s, when soil deposits become stable according to kinematic energy and unbalanced forces of the system (Kafaji 2013).

Table 1 Statistical properties of the soil parameters

Deterministic analysis of landslide runout motion

In this section, the example is modeled deterministically, i.e., as a homogeneous slope, in order to provide the benchmark for the subsequent analysis of the influence of spatial variability of soils on post-failure behavior. In this model, all mechanical parameters are constant; the mean values of $c$ and $\varphi$ are used as their peak values. Figure 5 shows the runout motions of the landslide computed with the deterministic model at critical times (colored contours represent the equivalent plastic strain invariant and the runout distance). At t = 3.0 s, the destruction of the failure block is evident; at this moment the mass slides about 1.7 m away from the initial toe of the slope (Fig. 5a). Also, the failure block is fully detached from the original slope and the rotational rupture surface is developed from the toe of the slope to the topside. At t = 5 s, the maximum soil displacement reaches 7.55 m (Fig. 5b). Then, the sliding mass further moves downward along a fully connected rotational rupture surface and spreads onto the ground. During the post-failure process, a large deformation can be observed, and the landslide is completely deposited at t = 11.0 s (Fig. 5c). For a quantification of the failure consequence, the runout distance with 10.78 m length is obtained from this calculation, which elaborates how long the failure mass moves from the initial toe of the slope to the new-formed deposit head.

Fig. 5

Configurations of landslide by deterministic analysis at critical times. a time at 3 s, b time at 5 s, c time at 11 s

Comparative analyses and stochastic assessment

In this section, comparative analyses are conducted in order to investigate the uncertainties of the post-failure behavior in the spatially varying cohesive-frictional deposit. The case with ${COV}_{c}$ = 0.1, ${COV}_{\varphi }$= 0.1, ${\rho }_{c, \varphi }$= − 0.5 is taken as a basic case in the following stochastic analysis.

Figure 6a and b show typical RF samples of soils’ $c$ and $\varphi$ values in the slope, respectively, which is randomly selected from samples. The blue color represents relatively weak parts with low cohesion or friction angle values, while the red color represents relatively strong sections with high cohesion or friction angle values. Given the negatively correlated coefficients, the increase of one parameter value decreases the other parameter’s value; for instance, the soil with a large cohesion would have a low friction angle. By conducting the simulation, the runout distance of the landslide in the heterogeneous slope with the spatially varying properties is obtained as 12.14 m (see Fig. 6c). Compared with the homogeneous analysis, the heterogeneous case resulted in a different runout distance. It should be noted that the results from different RF samples may differ due to the random distribution and combination of $c$ and $\varphi$. Therefore, multiple MC simulations are required to capture the range of possible post-failure deformations of the landslide. It is noteworthy that a particular distribution of a bivariate RF may result in stable slopes that do not lead to a failure/post-failure scenario. In this work, the analysis of these limited cases is overlooked as they are out of the scope of the study. Given that the main focus of the study is on analyzing post-failure response, only cases where a material points’ displacement larger than 1 m occurs are considered as sliding/failure cases.

Fig. 6

A typical sample of bivariate RF of cohesion and friction angle, with cross-correlation coefficient of − 0.5. a variations in $c$, b variations in $\varphi$, c corresponding calculated runout distance

In the following MCS analysis, the basic case is used to check the convergence criterion at first. Figure 7a shows the variability results of runout distance calculation for 1000 realizations. The variation of the runout distance and corresponding COV are plotted as functions of the number of MCS. The figure shows that 1000 MC simulations can produce reasonably stable results and give reliable estimates of runout distances. Figure 7b shows the time-runout distance curves of landslides by MCS through 1000 samples, which presents the calculated runout distance of landslides at different times from all samples. Unlike the deterministic analysis, the heterogeneous cases show dramatic uncertainties of runout motions during the sliding time, in which the maximum runout distance is 12.23 m and the minimum runout distance is 9.40 m. As for larger runout distance, this is mainly because in some cases the failure paths in the spatially varying soil generally go through relatively weak soil zones (Liu et al. 2019; Huang et al. 2020), which may lead to an extensive runout distance. While other cases may face strong soil sections through their failure paths and consequently have relatively smaller runout distances. Moreover, the runout distance of stochastic cases varies substantially at different times, and they distinctly differ from that of the deterministic analysis. For the same time, after 1.0 s, the runout distances of different samples start to differ from one case to another and vary in a relatively wide range of about 3 m at t = 3 s. It can be seen that for some cases the landslide does not initiate at this moment (i.e., runout distance = 0 m), while at the same time the runout distances of other cases can reach up to about 3 m. Although the runout distance varies from one realization to another, the landslides reach their final equilibrium state within a certain range of time from 6.0 to 8.0 s. Based on Fig. 7c, the variability of runout distance among all samples firstly increases with time up to t = 4.31 s. At this moment, the standard deviation of runout distance reaches to its peak value with 0.71 m, which means the maximum variance of runout distance can be observed at this time (see Fig. 7c). After this, the variance of runout distance decreases until it meets a stable state.

Fig. 7

a The convergence of calculated mean value and COV of runout distances; b time-runout distance curves for MCS samples and comparison with the deterministic analysis; c standard deviation of runout distances from MCS samples; d probability density histogram for the calculated runout distance fitted by a normal distribution

In order to reasonably estimate the distribution of computed runout distances, three common candidate distributions including Normal, Gumbel and Weibull distributions are examined. A goodness-of-fit test method is applied to identify the best-fit marginal distribution underlying computed data. Both the Akaike information criterion (AIC) and Bayesian information criterion (BIC) (Akaike 1974; Wang et al. 2020a) are adopted to identify the best-fit distributions; they can be denoted as:

$$\begin{array}{c}AIC=-2\sum\limits_{i=1}^{N}{\text{ln}}f\left[\left({x}_{i};p,q\right)\right]+2{k}_{1}\end{array}$$

(20)

$$\begin{array}{c}BIC=-2\sum\limits_{i=1}^{N}{\text{ln}}f\left[\left({x}_{i};p,q\right)\right]+{k}_{1}{\text{l}}{\text{n}}N\end{array}$$

(21)

where the ${x}_{i}$ ($i=\mathrm{1,2},...,N$) is the output runout distance; $N$ is the number of samples; $f({x}_{i};p,q)$ is the probability density function (PDF) of alternative marginal distribution function, and $p$ and $q$ are the parameters of distribution; ${k}_{1}$ is the number of distribution parameter in the alternative marginal distribution. For these abovementioned marginal distribution functions, the ${k}_{1}=2$ (Burnham and Anderson 2004). According to Eqs. (20) and (21), the corresponding AIC values for the Normal distribution, Gumbel and Weibull distribution are 1.52 $\times {10}^{3}$, 1.62 $\times {10}^{3}$, and 1.74 $\times {10}^{3}$, respectively; and their BIC values are 1.53 $\times {10}^{3}$, 1.63 $\times {10}^{3}$, and 1.75 $\times {10}^{3}$, respectively. Note that the distribution associated with the smallest AIC and BIC values is identified to be the best-fit marginal distribution to the output runout distance. To this end, both the AIC and BIC values indicate that the Normal distribution is the best-fit distribution for the output runout distance. Subsequently, the maximum likelihood method is adopted to estimate the corresponding mean value and standard deviation of the Normal distribution.

Figure 7d shows the probability density histogram for the calculated runout distance from the basic case with 1000 samples. In this figure, the horizontal coordinate represents the calculated runout distance, the left vertical coordinate represents the PDF, and the right coordinate represents the CDF. The red line is the PDF result fitted by the Normal distribution based on the computed values, and the blue line is the corresponding CDF result. It can be seen that the Normal distribution can relatively describe the predicted runout distances, with an estimated mean value of 10.80 m and standard deviation of 0.51 m. If the deterministic result is considered a limit of safety, the exceedance probability (runout distance > 10.78 m) is about 52%. Meanwhile, the largest runout distance with 12.23 m obtained from the stochastic analysis notably outstrip the limit of safety. Both results show a considerable discrepancy in post-failure motions, in terms of runout distance, between the homogeneous analysis and the heterogenous analysis. Consequently, consideration of the spatial variability of a cohesive and frictional soil significantly increased the uncertainty in estimation of the runout distances; that underpins the need to incorporate the effect of natural heterogeneity of soil in modeling of the post-failure processes.

Parametric analysis on the statistical information of soil properties

In this section, a parametric analysis considering different COVs of $c$ and $\varphi$ and their cross-correlation coefficients is performed to study the effect of the statistical information of soil properties on the results of probabilistic runout analysis. The cases with different parametric settings or auto-correlation structures have been checked by the same convergence criterion, and for each case it has been observed that 1000 MC simulations satisfies the convergence criterion.

Figure 8 shows the PDFs, CDFs, mean values, and exceedance probability of runout distances with five different ${\mathrm{COV}}_{c}$ ranging from 0.1 to 0.5. According to the results, the computed runout distance varies over a wide range from less than 8 m to over 15 m for all cases, which indicates that there is a large uncertainty in the landslide’s runout motions. It can be found that the shape of the PDFs is strongly affected by ${\mathrm{COV}}_{c}$, as the PDF curves become wider and have lower peaks with increasing ${\mathrm{COV}}_{c}$ (Fig. 8a). It indicates that higher variability of cohesion increases the probability of extreme cases (i.e., landslides with excessively long runouts). The variability of runout distance also increases with the increasing of ${\mathrm{COV}}_{c}$, and only little increments of corresponding mean values are observed in these cases. Figure 8d) shows the exceedance probability (${P}_{E}$) of runout distance for all cases, where horizontal coordinate represents the runout distance and vertical coordinate represents the ${P}_{E}$ that the runout distance is larger than a given distance. It is shown that in all cases the values of ${P}_{E}$ decrease as the runout distances increase. Among the ${\mathrm{COV}}_{c}$ cases, the influence of the ${\mathrm{COV}}_{c}$ on ${P}_{E}$ is significant. For example, if the acceptable risk probability is 0.1% (Wang et al. 2019), the corresponding runout distance of landslides with ${\mathrm{COV}}_{c}$ ranging from 0.1 to 0.5 are 12.24 m, 13.03 m, 13.5 m, 14.24 m, and 14.78 m, respectively. Therefore, any usage of lands within 14.78 m from the toe of the slope should be evaluated with caution if considering the worst case with large cohesive anisotropy.

Fig. 8

The effect of ${\mathrm{COV}}_{c}$ on PDFs, CDFs, mean values, and exceedance probability of runout distance. a PDFs, b CDFs, c mean values, d exceedance probability

Comparing with the features of the ${\mathrm{COV}}_{c}$ cases, the ${\mathrm{COV}}_{\varphi }$ in the slope not only influences the shape of the PDFs but also strongly affects the mean values of the runout distance (as shown in Fig. 9). It can be observed that with the increase of ${\mathrm{COV}}_{\varphi }$, the associated PDF curves become wider, lower, and farther to the left (Fig. 9a). The mean value and variance of the runout distance increase notably with the increase of the ${\mathrm{COV}}_{\varphi }$, which is consistent with the recent random SPH studies conducted by Zhang et al. (2020). For a small COV ${}_{\varphi }$, the varying range of $\varphi$ is correspondingly narrowed, and consequently the distribution of the runout distances will become narrower. Conversely, a large COV ${}_{\varphi }$ causes a wide range of values for $\varphi$, and subsequently results in a wider distribution of runout distances. Similar to the ${\mathrm{COV}}_{c}$, the influence of ${\mathrm{COV}}_{\varphi }$ on ${P}_{E}$ of runout distance is non-negligible. For example, if there is a transportation road located within a distance of 15 m away from the toe of the slope, the probability that it is affected by the landslide is 1.4% for ${\mathrm{COV}}_{\varphi }=0.5$. As for other ${\mathrm{COV}}_{\varphi }$ values, there is no possibility that the landslide would affect the road. Therefore, much larger runout distance is revealed by considering a large degree of heterogeneity in the friction angle property of the soil (i.e., a large value of ${\mathrm{COV}}_{\varphi }$).

Fig. 9

The effect of ${\mathrm{COV}}_{\varphi }$ on PDFs, CDFs, mean values, and exceedance probability of runout distance. a PDFs, b CDFs, c mean values, d exceedance probability

The effect of varying the cross-correlation between $c$ and $\varphi$ is also investigated. Figure 10 compares the PDFs, CDFs, mean values, and ${P}_{E}$ corresponding to different degrees of cross-correlation coefficient ${\rho }_{c, \varphi }$. Overall, it is observed that the shape and evolution of the PDF and CDF are less sensitive to the ${\rho }_{c, \varphi }$ compared with the ${\mathrm{COV}}_{c}$ and COV ${}_{\varphi }$ cases (Fig. 10a). In addition, the results illustrate that there are little differences among mean values. According to the figure, taller and narrower PDF curves are obtained when decreasing the ${\rho }_{c, \varphi }$ from 0 to − 0.7, which implies that the variation of the runout distance decreases with the decreasing of cross-correlation between $c$ and $\varphi$. As for negative coefficients, this could be partly explained by the negative correlation between $c$ and $\varphi$; the low values of cohesion are associated with high values of friction angle and vice versa. In other words, the uncertainty of soil shear strength estimation lowers when taking a larger negative correlation between $c$ and $\varphi$. Therefore, the variation of the total shear strength in soil is reduced and consequently the variance of runout distance is also reduced. As for cross-correlation equal to zero, the uncorrelated case shows biggest variance on runout distance due to bigger uncertainty caused by the uncorrelatedness between shear strength parameters. With regard to the ${P}_{E}$ (Fig. 10d), as the runout distance increases from 9 to 11 m, and the differences in ${P}_{E}$ curves related to different ${\rho }_{c, \varphi }$ cases are unnoticeable; however, the differences become apparent as the runout distance further increases (11 m to 14 m).

Fig. 10

The effect of ${\rho }_{c, \varphi }$ on PDFs, CDFs, mean values, and exceedance probability of runout distance. a PDFs, b CDFs, c mean values, d exceedance probability

Based on the above analysis, it can be concluded that the spatial variability and the negative cross-correlation of $c$ and $\varphi$ notably influence the post-failure behavior of landslides. A slope with larger ${\mathrm{COV}}_{c}$ would have a higher uncertainty of runout motions. While COV ${}_{\varphi }$ of soil in a slope have dual effect on both mean values and standard deviations of the runout distance. For the large COV ${}_{\varphi }$ case, it may result in longer runout motions and higher variance in the resulting runout distance. As for the cross-correlation between $c$ and $\varphi$, it has negligible effect on mean values of runout distance, but it has notable influence on their standard deviations, which increases the uncertainty of post-failure behaviors. Therefore, rational assessment of landslide post-failure behavior and its associated runout motions would need to be assisted with the appropriate parametric settings of the geotechnical bivariate RFs.

Effect of auto-correlation structures

Sensitivity analyses are performed on three auto-correlation structures, namely, isotropy, transverse anisotropy, and rotated anisotropy of soil shear strength parameters to investigate how anisotropically deposited soils influence the runout motions of a landslide, and to identify the critical auto-correlation structure which leads to the most extensive runout motion behavior of heterogeneous landslides.

Figure 11 presents typical realizations of the bivariate RFs of $c$ and $\varphi$ with different structures and their corresponding runout distances. For all cases, a commonly used cross-correlation coefficient is considered, as ${\rho }_{c, \varphi }$= − 0.5. The auto-correlation distances are set to ${\theta }_{1}$=20 m and ${\theta }_{2}$=2 m for all anisotropic cases. Figure 11a and show two isotropic cases (${\theta }_{1}$=${\theta }_{2}$= 2 m and ${\theta }_{1}$=${\theta }_{2}$= 20 m); and a typical realization of transverse anisotropy can refer to Fig. 6 in the previous section. Figure 11c-e consider rotated anisotropy, where the rotation angle $\beta$ is considered equal to − $45{ }^{\circ }$, $90{ }^{\circ }$, and $45{ }^{\circ }$, respectively. Obviously, for the rotated anisotropy cases, the slope is separated by inclined layers with different shear strength values, which has particular implications for the landslide failure mechanisms (Ma et al. 2018). According to Fig. 11, a considerable difference can be observed among different auto-correlation structures, in which the minimum runout distance is 8.99 m (Fig. 11c), while the maximum runout distance can reach up to 15.29 m (Fig. 11e). It should be noted that each realization of the generated RFs may be different, and therefore 1000 MC simulations for each case are conducted to reveal the post-failure behaviors of landslides for different auto-correlation structures. The mean value and variance of runout distances are obtained from the results of the MC simulations.

Fig. 11

Typical realizations of random fields of cohesion and friction angle and the corresponding runout distance. a isotropy (${\theta }_{h}$=${\theta }_{v}$=2 m), b isotropy (${\theta }_{h}$=${\theta }_{v}$=20 m), c rotated anisotropy ($\beta ={-45}^{^\circ }$), d rotated anisotropy ($\beta ={90}^{^\circ }$), e rotated anisotropy ($\beta ={45}^{^\circ }$)

To investigate the influence of the anisotropy and fabric orientation of soil shear strength parameters, 12 cases of anisotropy with rotated angles $\beta$ ranging from ${-75}^{^\circ }$ to ${90}^{^\circ }$ are considered, in which the $\beta ={0}^{^\circ }$ represents the transverse anisotropy case. According to all cases in Fig. 11, the runout distances of different landslide models may significantly differ when comparing any two different stochastic cases. However, after an average of a thousand MC simulation, the differences are gradually diminished. Figure 12 shows the variations of mean and standard deviation of runout distance of landslides as a function of the rotated angle $\beta$, in which two isotropic cases are plotted as a comparison. In this figure, the means and standard deviations of runout distance of the two isotropic cases are positively correlated to ${\theta }_{h}$, which is consistent with Zhang et al. (2020). This is because sliding mass is attracted to a pocket of weak soils and a large value of ${\theta }_{h}$ usually indicates a larger extent of a weak zone; it may result in larger deformation and failure behavior of landslides. Meanwhile, uncertainties of the corresponding response of landslides after failure would increase under the condition of high value of ${\theta }_{h}$, which is manifested as relatively high value of standard deviation. According to the figure, the estimated first and second moments of the two isotropic cases could preserve as upper bounds and lower bounds for rotated anisotropic cases (blue dotted lines in Fig. 12). Besides, the transverse anisotropy cases can serve as a boundary value for separating the rotated anisotropy cases.

Fig. 12

Influence of the auto-correlation structure on the runout distance of landslides. a mean values, b standard deviations

As for the rotated anisotropy cases, the mean values of runout distance for $\beta >{60}^{^\circ }$ and ${-75}^{^\circ }<\beta <{0}^{^\circ }$(equal to ${60}^{^\circ }<\beta <{180}^{^\circ }$) are less than the value of the transverse anisotropy case (black line) with 10.80 m ($\beta ={0}^{^\circ }$). This is the same for the standard deviations in these cases which are lower than the transverse anisotropy case with about 0.51 m. According to their estimates, these slope structures with relatively high angles for dipping layers or anti-dipping layers tend to form a smaller variability and short runout distance. When ${-75}^{^\circ }<\beta <{-30}^{^\circ }$, the mean values and standard deviations do not vary considerably but remain relatively constant among these cases. In addition, it can also be observed that the minimum mean value of the runout distance from the rotated anisotropy cases is close to 10.44 m which occurs at approximately $\beta ={-30}^{^\circ }$ and almost reaches to the result from the isotropic case (${\theta }_{1}$=${\theta }_{2}=$ 2 m) with 10.40 m. When the fabric orientation of soil layer is in the range of ${0}^{^\circ }<\beta <{60}^{^\circ }$, the corresponding mean values of runout distance are relatively higher and associated with more dangerous consequences of sliding, comparing with the transverse anisotropy case. Meanwhile, it can be concluded that the uncertainty of runout distance increases with the rotation angle from ${-30}^{^\circ }<\beta <{30}^{^\circ }$. The maximum mean value of runout distance with 11.35 m is obtained at $\beta ={30}^{^\circ }$, which is slightly lower than the isotropic case with 11.48 m (${\theta }_{1}$=${\theta }_{2}$=20 m).

According to Fig. 12, with respect to the mean and standard deviation values, there are two turning points at (a) $\beta ={30}^{^\circ }$ and (b) $\beta ={60}^{^\circ }$. The figure shows that that the fabric orientation of $\beta ={30}^{^\circ }$ in slopes is a critical angle that leads to the highest mean and standard deviation values for the runout distance. Furthermore, the standard deviations tend to accelerate increasing when the rotational angle is close to ${30}^{^\circ }$ which makes the uncertainty of runout motion higher. On the other hand, the figure also shows that at a limiting value of $\beta ={60}^{^\circ }$ the runout distances (both in terms of mean and standard deviation) switch from being higher than those relevant to transverse anisotropy case to comparatively lower amounts. The mean value and standard deviation for ${0}^{^\circ }<\beta <{60}^{^\circ }$ are higher than those for $\beta ={0}^{^\circ }$ (which represents a transverse anisotropy case, i.e., layered deposits). For $\beta >{60}^{^\circ }$ the mean value and standard deviation are smaller than those for $\beta ={0}^{^\circ }$. This finding is generally consistent with those from other relevant studies (e.g., Cheng et al. 2018; Huang et al. 2021) which estimated failure scale of slopes considering different rotation angles. Therefore, it becomes clear that fabric anisotropy is a crucial factor that can significantly and quantitatively influence the post-failure response of landslides.

Conclusions

In this paper, the effects of deposition anisotropy in cohesive-frictional soils on the potential landslide and post-failure behavior are analyzed using large strain GIMP simulation under an MCS framework with cross-correlated bivariate RFs for the strength parameters (i.e., cohesion and friction angle). Parametric analyses considering different COVs of $c$ and $\varphi$ and their cross-correlation coefficients are performed to study the effect of the statistical information of soil properties on the probabilistic runout behavior. The influence of orientation of anisotropic deposition on the runout motion of landslides is also investigated. Based on the results, the following conclusions can be drawn:

• 1. The spatial variability and negative cross-correlation of $c$ and $\varphi$ notably influence the post-failure behavior of landslides. Hence, for a realistic assessment of landslide post-failure behavior and associated runout motions, appropriate consideration of cross-correlation of $c$ and $\varphi$, needs to be taken into account.

• 2. Different auto-correlation structures of soil profile can result in different runout distances. Hence, it is essential to identify the appropriate auto-correlation structure of the deposition in assessing the large-deformation behavior of a landslide.

• 3. Rotation of soil deposition orientation significantly affects the runout motion. The deposition orientation at $\beta ={30}^{^\circ }$ in the analysis case is identified to produce the highest mean value and standard deviation of runout distance. The finding from this research highlights the significance of involving the orientation of soil stratification rather than the magnitude of shear strength alone in assessing the post-failure behavior of landslides.

Appendix. GIMP computation

In GIMP, a material body within the domain of the problem is discretized into a finite number of material points. The material point is a domain which mass and other variables vary according to a particle characteristic function ${\chi }_{p}(x)$ (Bardenhagen and Kober 2004) satisfying the partition of unity property ($\sum_{p}{\chi }_{p}(x)=1,\forall x$). In the particle domain ${\Omega }_{p}$, the particle characteristic function is generally been assumed

$${\chi}_{p}\left({x}\right)=\begin{cases}1\quad x\in {\Omega}_{p} \\0\quad otherwise\end{cases}$$

(22)

The function ${\chi }_{p}({\varvec{x}})$ is used as a basis for representing particle data by giving a material point property ${f}_{p}$ (e.g., particle density ${\rho }_{p}$, stress ${\sigma }_{ijp}$ and volume ${V}_{p}$).

To discretize the space occupied by the grid, computational grid shape function ${S}_{I}$ is introduced and required to be a partition of unity ($\sum_{I}{S}_{I}(x)=1,\forall x$), which is expressed as

$${S}_{Ip}=\frac{1}{{V}_{p}}\underset{{\Omega }_{p}\cap \Omega }{\int }{\chi }_{p}(x){N}_{I}\text{d}\Omega$$

(23)

$${S}_{Ip,j}=\frac{1}{{V}_{p}}\underset{{\Omega }_{p}\cap \Omega }{\int }{\chi }_{p}(x){N}_{I,j}\text{d}\Omega$$

(24)

where ${N}_{I}$ is the grid nodal shape function, and gradients ${S}_{Ip,j}$ are implicit functions of grid node position ${x}_{I}$, particle position ${x}_{p}$, and current particle volume $\Omega_{p}$. For this step, the discretization is analogous to FEM. After GIMP discretization, computing the total nodal force ${f}_{iI}$ is carried out

$$\begin{array}{c}{f}_{iI}={f}_{iI}^{int}+{f}_{iI}^{ext}\end{array}$$

(25)

$$\begin{array}{c}{f}_{iI}^{int}=-\sum_{p}{\sigma }_{ijp}{S}_{Ip,j}{V}_{p}\end{array}$$

(26)

$$\begin{array}{c}{f}_{iI}^{ext}=\sum_{p}{m}_{p}{S}_{Ip}{b}_{ip}+{\int }_{{\Gamma }_{t}}{S}_{I}\left(x\right){\overline{t}}_{i}\text{d}{\Gamma }_{t}\end{array}$$

(27)

where, the ${f}_{iI}^{int}$ is the nodal internal force and ${f}_{iI}^{ext}$ is the nodal external force at node I, respectively. Then, by integrating momentum equations and applying the boundary conditions, the material point position and velocities can be computed.