Rainfall-Induced Landslides from Initialization to Post-Failure Flows: Stochastic Analysis with Machine Learning

Abstract

Rainfall-induced landslides represent a severe hazard around the world due to their sudden occurrence, as well as their widespread influence and runout distance. Considering the spatial variability of soil, stochastic analysis is often conducted to give a probability description of the runout. However, rainfall-induced landslides are complex and time-consuming for brute-force Monte Carlo analyses. Therefore, new methods are required to improve the efficiency of stochastic analysis. This paper presents a framework to investigate the influence and runout distance of rainfall-induced landslides with a two-step simulation approach. The complete process, from the initialization of instability to the post-failure flow, is simulated. The rainfall infiltration process and initialization of instability are first solved with a coupled hydro-mechanical finite element model. The post-failure flow is simulated using the coupled Eulerian–Lagrangian method, wherein the soil can flow freely in fixed Eulerian meshes. An equivalent-strength method is used to connect two steps by considering the effective stress of unsaturated soil. A rigorous method has been developed to accurately quantify the influence and runout distance via Eulerian analyses. Several simulations have been produced, using three-dimensional analyses to study the shapes of slopes and using stochastic analysis to consider uncertainty and the spatial variability of soils. It was found that a two-dimensional analysis assuming plain strain is generally conservative and safe in design, but care must be taken to interpret 2D results when the slope is convex in the longitudinal direction. The uncertainty and spatial variability of soils can lead to the statistic of influence and runout distance. The framework of using machine-learning models as surrogate models is effective in stochastic analysis of this problem and can greatly reduce computational effort.

1. Introduction

Rainfall-induced landslides occur frequently around the world; they are a threat to life and cause huge economic losses [1,2,3,4,5,6]. According to a previous study [7], of the 4862 investigated landslides from 2004 to 2016, 79% were triggered by rainfall and led to the deaths of 55,997 people. Rainfall-induced landslides usually happen so rapidly that no mitigation measures can be introduced after instability is initiated [8,9]. Therefore, active measures (e.g., vegetations, retaining structures, piles, mesh, geocells, etc.) must be erected beforehand and the influence and runout distance must be determined, taking into account the probable soil-loss problems [10]. In addition, effective disaster management and slope stability analysis must be performed [4,11,12].

Rainfall-induced shallow landslides happen in two stages, as pointed out by Cascini [13]: first, the development of a complete shear band in the soil, which is termed the failure stage, and second, rapid post-failure flow, i.e., the post-failure stage. The second stage is the result of the failure stage. Therefore, the entire process of rainfall-induced landslides can be simulated according to two interrelated steps, which are rainfall infiltration analysis and post-failure large deformation analysis. The same two-step calculation has been validated with two case studies in China and Japan, respectively [14]. The fracture behavior that normally occurs in rocks is not considered in this study [15,16,17]. The finite element (FE) method, using coupled hydro-mechanical models, is usually adopted for slope stability analysis under rainfall infiltration [18,19,20] and is also employed in this study. Several numerical models have been developed to simulate the large deformation of post-failure flows, including but not limited to the particle finite element method [21], material point method [22], and smoothed-particle hydrodynamics [23]. The coupled Eulerian–Lagrangian (CEL) method is a technique used to model the flow of Eulerian materials through a fixed mesh by tracking the Eulerian volume fraction (EVF); it is also suitable for simulating the post-failure large deformation of landslides. An equivalent strength method is used to connect the two steps (infiltration analysis and post-failure flow analysis), which considers the variations in soil material properties caused by rainfall infiltration.

Most previous rainfall-induced landslide simulations employ deterministic analysis [24,25,26]. However, due to the limited number of tests conducted for most projects, soil properties cannot be precisely determined in practice, so stochastic analysis is often conducted using Monte Carlo simulations. Additionally, soil properties often vary spatially, due to the sediment history, and soil parameters are regarded as random fields [27,28,29,30]. Therefore, this study presents a framework for conducting stochastic analysis of runout distance and the influence distance of landslides, considering the spatial variability of soils. Additionally, three-dimensional simulations are conducted to study how the slope shape can influence the stability and runout of landslides.

One challenge of Monte-Carlo-based stochastic analysis with spatial variability is the high sampling demand, i.e., a large number of samples; therefore, simulations are needed for a single analysis, which demand considerable computing resources and time. Some efforts have been made to reduce the number of simulations [31,32]. Meanwhile, surrogate models and regression models are also used to replace the time-consuming numerical simulations [33,34]. The stochastic analysis framework was proposed as a result of training machine-learning models as surrogate models, using the framework to study slope stability [35].

The two-step calculations of rainfall-induced landslides are time-consuming in every simulation. However, an accurate probability density function (PDF) requires large sample sizes, which may require years to complete brute-force Monte Carlo analyses (directly simulated from the two-step calculations). The aim of this paper is to build a framework for a rainfall-induced landslide problem to significantly improve the efficiency of stochastic analysis: for this purpose, a small number of simulations are conducted to obtain the influence and runout distance. Then, the random fields of soil parameters and the calculated influence and runout distance are treated as the input and output to train a machine-learning (ML) model. This model will be used to predict the influence and runout distance for many samples and, thus, to estimate the PDF.

The structure of this paper is as follows. Soil models, the hydro-mechanical model, the CEL method, and, particularly, the method used to connect the infiltration analysis and post-failure flows are explained in Section 2. A deterministic analysis and parametric analysis are presented in Section 3, while three-dimensional analyses are conducted to examine how three-dimensional shape effects the accuracy of predictions. A stochastic analysis is presented in Section 4, including the generation of random fields, brute-force stochastic analysis, a brief introduction to neural networks, and stochastic analysis with machine learning.

2. Method

In this study, the rainfall infiltration process and the initialization of instability are simulated with a coupled hydro-mechanical FE model. The coupled Eulerian–Lagrangian method is used for post-failure flows.

2.1. Soil Constitutive Model

The Mohr–Coulomb model has been widely used in geotechnical engineering because (i) it has a small number of parameters that can be determined easily, and (ii) the concept is simple and can reflect the characteristics of both frictional and cohesive materials, such as soils. However, its yield surface does not employ continuous derivatives, which leads to the difficulty of convergence in some simulations. The Drucker–Prager model can, then, be used as an alternative way to avoid this problem, especially in large deformation simulations [36,37]. In this study, an extended Drucker-Prager model with a non-associated flow rule is used in both the first step (rainfall infiltration) and the second step (post-failure flows), and the yield function is:

$$f=t+p \tan \beta -d$$

(1)

$$t=\frac{q}{2} \left[1+\frac{1}{k_t}-\left(1-\frac{1}{k_t}\right){\left(\frac{J_3}{q}\right)}^3\right]$$

(2)

where $t$ is the effective shear stress; $J_2$ and $J_3$ are the second and third invariants of deviatoric stress; $q$ is generalized shear stress; $p$ is the mean stress; $k_t$ is the ratio of the yield stress in triaxial tension to that in triaxial compression and it controls the shape of the yield stress in the $\pi$ plane. When $k_t=1$ (as used in this study), this reduces to the convectional Drucker–Prager model where $t=q=\sqrt{3J_2}$. The Drucker–Prager friction angle, $\beta ,$ is the slope of the linear yield surface in the $p$–$t$ plane and is related to the true friction angle of the material; the Drucker–Prager intercept, $d,$ is related to the true cohesion and friction angle of the material. Figure 1a shows the linear Drucker–Prager model in the meridional plane, while Figure 1b shows a comparison of the linear Drucker–Prager model ($k_t=1$) with the Mohr–Coulomb model ($\phi =20°$).

The flow potential, $G,$ is:

$$G=t-p \tan \psi$$

(3)

where $\psi$ is the dilation angle.

2.2. Coupled Hydro-Mechanical Analysis

The soil is a three-phase mixture with the soil skeleton, water, and air. The total stress ${\sigma }_t$ is related to the water pressure $u_w$, air pressure, $u_a,$ and the effective stress, ${\sigma }^{\prime }$. If the air inside and outside is connected, the air pressure can be ignored, and the effective stress principle is:

$${\sigma }^{\prime }={\sigma }_t-\chi u_w I$$

(4)

where $\chi$ is Bishop’s parameter, which is approximately the value of the saturation in previous studies [38,39]. According to Arifin and Schanz [38], $\chi$ (ranging from 0 to 1) can be measured in the laboratory and $\chi =1$ means full saturation; $I$ is the second-order identity tensor.

Darcy’s law is expressed as:

$$sn{v}_w=-k_d \frac{\partial h}{\partial x}$$

(5)

$$h=z+\frac{u_w}{\left|g\right| {\rho }_w}$$

(6)

$$k_d=k_s \overline{k_d}$$

(7)

where $s$ is the saturation; $n$ is the porosity; $v_w$ is the seepage velocity of water; $k_d$ is the permeability of the soil; $h$ is the hydraulic head in the soil; $x$ is the coordinate; $z$ is the elevation above the reference elevation; $g$ is the acceleration of gravity; ${\rho }_w$ is the density of the water; $k_s$ is the relative permeability and its cubic power of saturation for the uniform pore-size distribution [40]; $\overline{k_d}$ is the permeability of fully saturated soil.

Porosity has a significant effect on slope stability. Soils with greater porosity have a larger water-retention capacity; therefore, for a certain rainfall intensity, slopes take longer to become saturated and, thus, unstable [41]. This phenomenon can be reproduced with the present model. However, some effects are hard to implement in this kind of two-step analysis, and are, therefore, ignored. For example, the soil–water characteristic curve is greatly influenced by porosity [42], which is not considered in the FE model. Some researchers [43,44] have considered the spatial variability of the porosity and proved that this variability may lead to some unexpected effective stress distribution and may make the failure process more complicated. Our stochastic analysis considers the spatial variability of soil strength parameters and permeability, but not their porosity. In the post-failure stage, the porosity of soils will also undergo dramatic change because of the rearrangement of soil particles and segregation [45,46]. The variation of porosity is difficult to establish, especially considering the effect of segregation [14]. In this study, porosity is assumed to be constant in post-failure flows.

The Abaqus FE software [37] with a coupled hydro-mechanical model is used to simulate rainfall infiltration. With the rainfall continuing, the deformation of the slope may increase rapidly, leading to instability at some stage, and finally, to landslides. This post-failure flow is large deformation, in which a traditional FE simulation may suffer mesh distortion and, thus, non-convergence. The results of the coupled analysis after slope failure would be mapped into the CEL model as the initial conditions.

2.3. Coupled Eulerian–Lagrangian (CEL) Method

In the CEL method, a mesh is fixed for the Eulerian domain, and materials such as soils can move freely in the mesh, which can help to avoid the mesh distortion problem near slip surfaces (where the shear behavior concentrates). Each iteration of the CEL model is illustrated in Figure 2. The configuration of a specific material in the Eulerian domain is tracked by its Eulerian volume fraction (EVF). A value of EVF = 1 for an element means that this element is totally occupied by the materials. The sum of EVFs in each element cannot exceed 1, and EVF = 0 means that this element is empty. The first step in the CEL method is similar to an updated Lagrangian FE simulation (Figure 2a); in the second step of the CEL method, as shown in Figure 2b (sometimes termed the Eulerian step), the mesh is reset, and a transfer algorithm is used to update all the variables and EVFs (Figure 2c).

In the CEL model, the contact between Eulerian domains and Lagrangian domains is modeled using the general contact method, which is based on the penalty method. In this study, the boundaries are modeled as Lagrangian rigid bodies. Seeds are created on the Lagrangian element faces and edges, while anchor points are created on the Eulerian material surface. The penalty method approximates spring deformation. The contact force, $Fp$, which is enforced between seeds and anchor points, is related to the penetration distance, $dp$, as:

$$F_p=k_p d_p$$

(8)

where the factor $k_p$ is the penalty stiffness, which depends on the Lagrangian and Eulerian material properties.

2.4. Transition between Finite Element Analysis and Coupled Eulerian–Lagrangian Analysis

The two steps should have the same physical conditions; the last increment of the first step will be the initial conditions of the second step. Therefore, the results exported from the coupled hydro-mechanical FE analysis will be imported into the CEL simulations. The exported data should include the coordinates, saturation, void ratio, stresses, and strain of each node or element, which is implemented by Python scripts. The nodal displacement and velocity will be ignored because they are too small, compared with the geometry of the slope.

In infiltration analysis, data are saved on the element nodes or integration points, and the deformed mesh has irregular shapes so that the meshes between the two steps are inconsistent. Therefore, the data from the first step cannot be used in the second step directly. The biharmonic spline interpolation method [47] is adopted to remap the data and is implemented using MATLAB scripts:

$$\mathrm{W}\left(x_i\right)={\displaystyle \sum }_{j=1}^N{\alpha }_j {\varphi }_m \left(x_i-x_j\right)$$

(9)

where $x_i$ and $x_j$ are the target interpolation points and initial data points, respectively; $\mathrm{W}\left(x_i\right)$ is the data on each target interpolation point; ${\alpha }_j$ is found by solving the linear system with all the known element nodal data; ${\varphi }_m$ is the biharmonic green function for each dimension and can be checked against previous studies [47].

Next, we substitute the effective stress of unsaturated soil (Equation (4)) into the yield strength (Equation (1)):

$${\tau }_f=d+{\sigma }^{\prime } \tan \beta =d-\chi u_w \tan \beta +{\sigma }_t \tan \beta .$$

(10)

As reported in previous studies [48,49,50], landslides often occur rapidly, meaning that excess pore water pressure does not have time to dissipate. This situation is similar to undrained conditions in soil mechanics [51]. Therefore, after the infiltration analysis, the pore water pressure is assumed to keep constant in post-failure flows, and an equivalent strength method is used [14]. For the soil in each FE element, the equivalent friction angle ${\beta }_e$ and equivalent cohesion $d_e$ used in CEL are:

$$d_e=d-\chi u_w \tan \beta$$

(11)

$${\beta }_e=\beta$$

(12)

Here, the pore water pressure, $u_w,$ and $\chi$ (i.e., saturation) are from the last iteration of the FE analysis. Due to the increase in water pressure and saturation, the equivalent cohesion is reduced, compared with that before rainfall occurred.

The equivalent friction and cohesion are based on the total stress; this stress is related to the equivalent density, ${\rho }_e$, which is increased due to the rainfall:

$${\rho }_e={\rho }_s+n s {\rho }_w.$$

(13)

3. Rainfall-Induced Landslide and Post-Failure Flow

Figure 3a,b shows the geometry of two types of slopes, which have the same sizes but are different in terms of layers. First, the uniform slope is studied. The rise and run of the slope are 10 m and 14 m, respectively. The uniform rainfall boundary condition is enforced on the top surface, while a completely fixed boundary condition is applied to the bottom. The left and right boundaries are fixed in the normal direction.

Table 1. Parameters used in the simulation of rainfall-induced landslides.

Parameters	Values
Young’s modulus $E$	100 MPa
Poisson’s ratio $v$	0.3
Cohesion $d$ in Drucker–Prager	10 kPa
Friction angle $\beta$ in Drucker–Prager	35°
Dilation angle $\psi$	0°
Soil particle density ${\rho }_s$	2650 kg/m3
Water density ${\rho }_w$	1000 kg/m3
Initial porosity $n$	0.3
Rainfall intensity $q_r$	0.018 m/h
Duration	20 h
Initial matric suction	20 kPa
Hydraulic conductivity $k$	0.036 m/h
SWCC parameter ${\alpha }^{\prime }$	0.31 m−1
SWCC parameter $n^{\prime }$	1.19
Hydraulic conductivity parameter $\eta$	1.962 m−1

shows the parameters used in the simulation of the base model. The soil–water characteristic curve (SWCC) is fitted using the van Genuchten model [52]. The relative permeability is modeled using the Gardner model [22,53]. Further simulations may only vary one parameter at a time and may fix the others. Figure 4 shows the SWCC and hydraulic conductivity that are used in simulations.

The FE mesh contains 440 elements, and the element size is approximately $1 \mathrm{m}\times 1 \mathrm{m}$. Each element contains four Gaussian integration points. First, a steady-state step is conducted to obtain an initial stress distribution without rainfall. Then, transient analysis is conducted to simulate the rainfall-induced hydro-mechanical response of the slope.

The variations in saturation with time are presented in Figure 5. In total, 40 h are simulated (20 h of rainfall and 20 h after rainfall) to check the distribution of saturation. The initial saturation of the slope is uniformly 0.684. After the start of rainfall, the saturation increases rapidly near the surface of the slope, forming a high saturation band (Figure 5a). With time, the high saturation band widens (Figure 5a,b) due to the rainfall. After the rainfall stops, this band will move downward and spread out with a decrease in saturation (Figure 5c).

The predicted displacement from the first step is presented in Figure 6. Figure 6a shows the rainfall of 10.6 h, the slopes have almost no deformation. Figure 6b shows the initial deformation of the slope in 18 h and in Figure 6c the sliding surface can be clearly observed. A scale factor is taken as 100 to amplify the displacement of the slope. Because damping is used, the model can cope with a certain degree of deformation without a failure of convergence. Figure 7a shows the displacement of the slope top and Figure 7b shows that of the slope toe. Figure 6 and Figure 7 illustrate that this failure happens at around 17 h. If the rainfall duration is less than 17 h, the slope is expected to remain stable. In addition, the displacement in the first step is quite small compared with the geometry of the slope; therefore, ignoring this displacement in the second step is reasonable.

A finer mesh size of $0.5 \mathrm{m}\times 0.5 \mathrm{m}$ is also employed in the first step and the results are shown in Figure 8. Figure 8a,b shows the saturation distribution and displacement in the FE infiltration analysis. The displacement is also amplified by 100 to clearly show the deformation of the slope. According to Figure 8, there is no significant difference between the fine and coarse mesh sizes. Taking into account the computational efficiency of the simulations, a mesh size of $1 \mathrm{m}\times 1 \mathrm{m}$ is sufficient to address this problem.

Each element in the CEL method contains only one Gaussian integration point. The CEL simulation domain contains 2400 elements, and the mesh size is $0.5 \mathrm{m}\times 0.5 \mathrm{m}$, which is consistent with the integration points in the first step. Initially, 1502 elements contain materials. Additionally, because the pore water pressure and saturation are different at different locations after infiltration, the equivalent cohesion $d_e$ in the CEL method are spatial variables, even if the initial slope has uniform strength parameters; 1502 different materials are filled in the CEL, and each is tracked by its own EVF.

Figure 9 shows the deformation of the slope. The slope does not register the tiny amount of velocity in the first step (Figure 9a). An initial sliding surface is clearly indicated in Figure 9b. The soil slides downward along this sliding surface (Figure 9c). The sliding mass gradually reduces (Figure 9d) and finally stops (Figure 9e).

The determination of runout and influence distance is illustrated in Figure 10. The dots show the sum of all 1502 EVFs for all CEL materials. Elements that are fully occupied by materials have values of EVF = 1, and empty elements have EVF = 0. Therefore, the final slope profile exits at elements that have 0 < EVFs < 1. The Eulerian analysis cannot track the material interface exactly. Instead, it can only be approximately recovered by conducting curve fitting (Figure 10). In order to obtain an accurate influence and runout distance, the top profile and bottom profile are fitted separately (Figure 10). Additionally, because the landslide front is represented by only one or two layers of elements, it is even harder to accurately obtain the runout distance; the different influence distance $D_{*}^i$ and runout distance $D_{*}^r$ are defined at different heights (Figure 10) to depict the final slope profile. Here, the superscript i denotes the influence distance, and r denotes the runout distance. Subscript numbers (0, 0.5, 1, 1.5) indicate the vertical distance between the defined height and the initial top and toe of the slope (Figure 10). Figure 11 shows the variations in influence distances (Figure 11a) and runout distances (Figure 11b) over time, which clearly shows that the post-failure flow stops at around 5 s. After 5 s of deformation, the slope head and slope toe stop moving and are finally stable. Additionally, the influence distance is generally larger than the runout distance, which is to be expected because only part of the upper soil eventually accumulates at the toe of the slope.

3.1. Effect of Soil and Rainfall Parameters

Rainfall-induced landslides are affected by many factors, including soil parameters, initial saturation, rainfall intensity, rainfall duration, and slope geometry. Therefore, it is necessary to conduct a sensitivity analysis, which is achieved by simply varying one parameter but fixing the other parameters in this study. In this section, seven parameters are studied, including the Drucker–Prager cohesion, $d$, the Drucker–Prager friction angle, $\beta$, soil particle density, ${\rho }_s$, hydraulic conductivity, $k$, initial saturation, $S_i$, rainfall intensity, $q_r$ and rainfall duration, $T_r$.

The variation in final influence and runout distance with each different Drucker–Prager cohesion $d$ (from 5 to 100 kPa) in the uniform slopes is presented in Figure 12a,b. When the cohesion is less than 5 kPa, the slope is initially not stable and will fail at the geostatic step. When the cohesion increases (higher strength), the influence distance and runout decrease (i.e., the slope is more stable). When the cohesion is greater than 50 kPa, this rainfall intensity and duration will not cause instability and landslides.

Soil properties usually vary greatly in the vertical direction due to sedimentation history. In this section, layered slopes are also considered (three layers, as in Figure 3b; the parameters for the different layers are listed in

Table 2. Parameters used for the layered slopes.

	Test Label	Top	Middle	Bottom	Mean	Max
Cohesion $d$ (kPa)	D1	5	7.5	10	7.5	10
	D2	7.5	10	12.5	10	12.5
	D3	5	10	15	10	15
	D4	15	20	25	20	25
Friction angle $\beta$ (°)	F1	31	33	35	33	35
	F2	34	35	36	35	36
	F3	33	35	37	35	37
	F4	35	37	39	37	39
Soil particle density ${\rho }_s$ (kg/ m3)	R1	2250	2450	2650	2450	2650
	R2	2550	2650	2750	2650	2750
	R3	2450	2650	2850	2650	2850
	R4	2650	2850	3050	2850	3050

). For example, to study slopes with different cohesion values, $d,$ for different layers, four simulations are conducted. Due to the sedimentation history, the bottom soils usually have greater strength and density, which pattern is followed in the simulations, as in Table 2. Tests D2 and D3 have the same average cohesion as the base model, while D3 has greater variation. Test D1 has smaller average cohesion than the base model, and test D4 has greater average cohesion. Figure 12c,d gives the influence and runout distance for these layered slopes; the horizontal axis is the maximum cohesion for the tests and so, from left to right, tests D1 to D4 are shown. The horizontal dash lines represent the results for the uniform slope (the base model). In Figure 12c,d, the variation of cohesion between layers has limited influence (i.e., D2 and D3 are the same), but the layered slopes (D2 and D3) demonstrate smaller runout distance than the uniform slope (the base model) because of the greater cohesion at the slope toe in D2 and D3. Additionally, in agreement with the trend of increasing cohesion in uniform slopes, the increase in average cohesion in the layered slopes (from D1 to D4) also leads to increased stability and reduced runout.

The results concerning the different Drucker–Prager friction angles, $\beta$ (from 32 to 45°), in the uniform slopes are illustrated in Figure 13a,b. The slope will initially be unstable and will fail in the geostatic step if the friction angle is less than 32°. When the friction angle increases (higher strength), the influence and runout distance decrease (i.e., more stable), and a friction angle greater than 45° will not have landslides. Similarly, four-layered slopes are studied with different friction angles, $\beta ,$ for the different layers (parameters are chosen based on the same logic as the study of cohesion), and Figure 13c,d gives the results. Tests F2 and F3 have smaller runout values than the base model (these three tests have the same average friction angle) and F3 is smaller than F2, which is because the larger friction angle at the slope toes is effective in reducing the runout. Additionally, a greater average friction angle (from F1 to F4) will also reduce the runout.

The results concerning particle density, ${\rho }_s,$ in the uniform slopes are illustrated in Figure 14a,b (with values from 1450 to 2950 kg/m³). When the density increases, the influence and runout distance increase (i.e., they become unstable due to the higher gravity load). Their impact on influence distance is very small (Figure 14a) but the density has a great impact on the runout (Figure 14b). Table 2 lists four tests of layered slopes with different particle densities; the values are not far from the typical values for sands and clays (2650 kg/m³), the obtained runout and influence distance only vary slightly (Figure 14c,d), and this variation is possibly smaller than the error associated with the determination of runout.

The impact of hydraulic conductivity, $k$ (from 0.018 to 0.7 m/h), is presented in Figure 15. To ensure that no runoff happens on the slope surface, $k$ cannot be less than 0.018 m/h (i.e., the rainfall intensity, $q_r$). When the hydraulic conductivity increases near the rainfall intensity (0.018 m/h), the runout distance will exhibit a violent drop (Figure 15b). When the hydraulic conductivity continues to increase, the influence and runout distance decrease (i.e., they become more stable), as shown in Figure 15a,b. This is because higher conductivity will let water infiltrate quickly through the slope and reach the bottom, and the growth of saturation in the slope will be reduced.

Initial saturation, $S_i$, is also a factor that influences the runout. The results are presented in Figure 16a,b ($S_i$ is from 0.512 to 0.83). When the $S_i$ is greater than 0.83, the slope is initially not stable and fails in the geostatic step. When the $S_i$ increases, the influence and runout distance increase slightly.

The rainfall intensity (varying from 0.012 to 0.036 m/h) and duration (from 1 to 40 h) were investigated, and the results are presented in Figure 17a,b and Figure 18a,b. The intensity value must be smaller than the hydraulic conductivity value (0.036 m/h), otherwise, surface runoff will occur. Figure 17 shows that when the rainfall intensity increases, the influence and runout distance also increase. In particular, when the rainfall intensity is very close to the maximum possible intensity and is greater than 0.03 m/h, the runout distance will rise dramatically (Figure 17b), which is similar to the drop seen in Figure 14b; in both cases, the rainfall intensity is very close to the hydraulic conductivity. In terms of the impact of rainfall duration, $T_r$, Figure 18a shows the different influence distances with changing $T_r$ and Figure 18b gives the results of the different runout distances. According to Figure 18, the increase in $T_r$ leads to an increase in influence and runout distance. However, the influence distance is not very sensitive to this variable, while the runout distance increases gradually with persistent rainfall.

3.2. Effect of Slope Shapes

The previous studies mostly employ two-dimensional (2D) analysis by assuming that the longitudinal length, perpendicular to the slope cross-section, is long, and the problem is simplified into a plane-strain problem. In this section, three-dimensional (3D) slopes are studied, including two types of geometry: a concave slope and a convex slope (Figure 19a,c). The side views of the two slopes are the same as those of the base model. In both tests, the Eulerian domain is discretized into 24,000 elements. Of these, 19,039 elements are initially occupied by materials for the concave slope, and 11,441 elements are occupied by materials for the convex slope. Because tracking tens of thousands of materials and their EVFs requires a large amount of RAM and CPU resources, each initially occupied element will be not modeled as a separate material, as in the 2D simulations; however, every 12 elements (concave) or 8 elements (convex) initially share the same material and are tracked by an EVF. Therefore, 1749 (concave) and 1490 (convex) materials are defined in the CEL. The material parameters are calculated via the same biharmonic spline interpolation method used in the infiltration simulations. The top view of the final profiles is shown in Figure 19b,d. The dash lines represent the contours at the same height as in the initial profile, with magenta for the top contour, red for the mid-segment, and blue for the bottom. The solid lines use the same contours for the final profile, while the dotted lines are drawn from the 2D predictions (i.e., the base model). The black arrows indicate the sliding directions. In both cases (concave and convex), the runout distances are smaller than in the 2D prediction; in particular, the runout distances at the corner are almost negligible. Therefore, the runout distance predicted by the 2D analysis is conservative and leads to safe design in practice. In terms of the influence distance, for those sections away from the corner, the predictions from 3D and 2D analyses are almost identical, with a 1.5% difference. For the concave slope, the top corner does not move, so the conservative prediction from the 2D analysis is safe in practice. However, in the case of the convex slope, the influence distance at the top corner is larger than in the 2D analysis. Therefore, great care must be taken to interpret the 2D analysis when the slope is convex in the longitudinal direction.

4. Stochastic Analysis with Machine Learning

Soil parameters are often associated with uncertainty and spatial variability. The proposed two-step framework is used to conduct stochastic analysis, wherein the initial parameters, including cohesion, $d$, friction angle, $\beta$, soil particle density, ${\rho }_s,$ and hydraulic conductivity, $k,$ are modeled with random fields. The material parameters are assumed to follow log-normal distributions. The exponential autocorrelation function ($\rho \left(x,x^{\prime }\right)=\exp \left(-\frac{\left|x-x^{\prime }\right|}{l_H}-\frac{\left|x-x^{\prime }\right|}{l_v}\right)$) is used, where $l_H$ and $l_v$ are the horizontal and vertical correlation lengths. These correlation lengths reflect the rate at which the correlation is delayed between two points in space. In other words, soil particles will be more similar with a shorter distance from each other. In the infiltration analysis, there are 440 elements; therefore, 440 materials are defined, which are spatial variables. The open-source software GSTools [54] was used to generate random field samples. Figure 20 shows the random fields of cohesion in the simulations. The parameters used in this stochastic analysis are listed in

Table 3. The parameters for the stochastic analysis.

Parameter	$\mathbf{Mean} \mathit{\mu }$	COV	${\mathit{l}}_{\mathit{H}}$	${\mathit{l}}_{\mathit{v}}$
Cohesion $d$	10 kPa	0.1	16 m	8 m
Friction angle $\beta$	35°	0.05	16 m	8 m
Soil particle density ${\rho }_s$	2650 kg/m3	0.05	16 m	8 m
Hydraulic conductivity $k$	0.036 m/h	0.3	16 m	8 m

. At this point, 1000 random fields are generated, and the corresponding simulations are performed. Each two-step simulation in this study is conducted on a laptop computer, with an Intel Core 7 CPU @ 2.80GHz and 16 GB of random-access memory (RAM). Each simulation takes 20–30 min to complete, meaning that 1000 simulations cost 16 days of computation.

The numerically estimated PDF of the influence and runout distance are illustrated in Figure 21 and Figure 22, respectively. The mean values of influence distance are larger than the runout distance, which is consistent with the results of the deterministic analysis. In addition, the standard deviation of influence distance is smaller than that of runout distance, which means that the uncertainty and spatial variable of soil parameters have a stronger impact on the uncertainty of runout distance than the influence distance. These simulations took over two weeks and the estimated PDFs are not very accurate, compared with the solid lines shown in Figure 21 and Figure 22, which were obtained from machine-learning-aided stochastic analysis using ${10}^5$ Monte Carlo samples.

Machine learning (ML) algorithms can build mathematical models, based on existing sample input-output pairs [35]. For the rainfall-induced landslides presented in this study, the inputs are the material parameters, which are spatially variable, and the outputs are the influence and runout distance. In the framework of machine-learning-aided stochastic analysis, the input data (spatially variable material parameters) and calculated influence and runout distance from a small number of two-step simulations are fed into the ML algorithms as training data. A general mathematical relationship will then be found and can be used as a surrogate model to predict the influence and runout distance, which is more effective than the two-step simulations.

An artificial neural network (ANN) is used in this study, which is inspired by the biological neural networks that constitute animal brains [35]. It is similar to the human brain’s neural network, from the perspective of information processing. A neural network is an operation model, which is composed of a large number of nodes (or neurons) connected with each other. These artificial neurons receive the signal, process it, and pass the signal on to adjacent neurons. ANNs can have multiple input and output connections and the connections between neurons are called “edges”. Signals are transferred from the input layer to the output layer. The connection between every two nodes represents a weighted value for the signal passing through the connection, which is called weight and is similar to the memory of the artificial neural network. Each node has a specific output function called the activation function. Before the output is produced, the weighted sum is combined with this neuron’s internal state (or activation) by using this activation function. Non-linear activation functions can help ANNs to learn more complex data and give exact predictions. The loss function is used to define the accuracy and ANNs are trained to minimize the loss function. The learning rate defines the size of the corrective steps that the model must take to adjust for errors in each observation. A high learning rate shortens the training time but may cause an unstable training process and may lead to a local minimum other than a global minimum, while a lower learning rate results in a long training process. The open-source library, TensorFlow, which implements the ANNs is used in the present study.

The FE model in infiltration analysis contains 440 elements and 4 random fields. Therefore, the input size of ML algorithms is 1760 and the output size is 8, which are the influence and runout distance defined at different heights. For this analysis, 1000 input-output pairs were obtained from two-step simulations, of which 500 samples were first used as training samples and 500 were validation samples. It was shown that deep learning neural networks can perform very complex tasks [55], but in the case of the simple problem in the present study, ANNs with only one hidden layer were enough. The activation function is the rectified linear unit (ReLU), while the mean absolute percentage error (MAPE) was chosen as the loss function. The Adam learning-rate optimization algorithm [56] was used and a learning rate of 0.001 was adopted in this study. An early-stop technique (with a maximum of 5000 epochs of training) was also used to stop the training after the loss function did not change.

The number of neurons is a tuneable hyperparameter in ANNs. Figure 23 shows the variety of MAPEs possible when the neurons change from 20 to 100. The increase in neurons will lead to more complex models and is expected to improve the accuracy. However, when the neurons are more than 60, no significant improvement can be observed. Therefore, an ANN with 60 hidden neurons is optimal. This size of ANN can handle very complex problems because it has more than 10⁵ trainable parameters. Figure 24 and Figure 25 compare the influence and runout distance calculated with two-step simulations with those predicted from machine-learning surrogate models. It is clear that the mean absolute percentage error (MAPEs) of runout distance (about 9%–10%) is greater than that of influence distance (4%–8%). This is because in Eulerian simulations, the interface and the profile cannot be exactly tracked but are instead recovered from the positions of elements where the 0 < EVFs < 1, using curve fitting. Therefore, there are errors when determining the runout or influence distance and this noise in the input data prevents the ML algorithms from producing very accurate models. The errors here (of around 10% for the runout distance and 6% for the influence distance) from machine-learning surrogate models are, thus, considered acceptable.

With this ANN as a surrogate model, the influence and runout distance for the millions of Monte Carlo samples can be easily evaluated. In Figure 21 and Figure 22, the black solid lines are the PDFs of influence and runout distances for the 10⁵ Monte Carlo samples, which are evaluated with the machine-learning surrogate model. The computation time for 500 simulations is about 7–10 days, while the time needed for the training and prediction of 10⁵ samples using ML is only about 10 min.

Table 4. The statistics of influence distance and runout.

Random Variable of Interest	Brute-Force Analysis		Machine-Learning-Aided
	Mean	Standard Deviation	Mean	Standard Deviation
$D_0^i$	7.72 m	0.72 m	7.85 m	0.66 m
$D_{0.5}^i$	5.40 m	0.60 m	5.32 m	0.55 m
$D_{1.5}^i$	3.81 m	0.58 m	3.86 m	0.49 m
$D_0^r$	4.79 m	1.02 m	4.88 m	0.97 m
$D_{0.5}^r$	3.80 m	0.82 m	3.82 m	0.80 m
$D_1^r$	3.32 m	0.73 m	3.37 m	0.68 m
$D_{1.5}^r$	3.00 m	0.65 m	3.07 m	0.65 m

shows the mean and standard deviations for the influence and runout distance. The errors between the brute-force stochastic analysis and machine-learning-aided stochastic analysis are smaller than 2%. However, machine-learning-aided analysis is able to obtain an accurate PDF, due to the larger number of samples used. Therefore, this framework greatly reduces the amount of computation necessary and ensures higher accuracy with a large number of Monte Carlo samples.

After the PDFs of influence and runout distance are obtained, the probability that an infrastructure (i.e., a house, railway, etc.) lies within the influence of landslides can be estimated using the following function:

$$p_f={{\displaystyle \int }}_a^{\infty }{f}_{PDF}\left(x\right) dx$$

(14)

where $a$ is the distance between the infrastructure and the slope top and toe, and $f_{PDF}\left(x\right)$ is the probability density function of the influence and runout distance.

5. Conclusions

This paper presents a study of the influence and runout distance of rainfall-induced landslides. In particular, the effect of the 3D slope shape was studied and stochastic analysis was conducted.

A two-step approach was adopted for rainfall-induced landslides. This approach can simulate the entire process, from initialization to post-failure flows. A coupled hydro-mechanical FE model was used to simulate the infiltration of rainfall, which can evaluate the pore pressure, displacement, and plastic strain of the slope. The coupled Eulerian–Lagrangian (CEL) method was used as the second step to simulate the post-failure flows. An equivalent strength method was used to build a connection between the infiltration analysis and the post-failure flow. A method was also developed to more accurately quantify the influence and runout distance using Eulerian analysis.

The sensitivity analysis showed that the influence and runout distance are affected by many factors, including the soil strength parameters (cohesion and friction angle), soil permeability and density, initial slope saturation, and rainfall intensity and duration. Two 3D slopes were studied with different shapes in the direction perpendicular to the cross-section (a concave slope and a convex one); it was found that the runout predicted from the 3D analysis is smaller than in the predictions of the 2D analysis, assuming plain strain; therefore, 2D analysis is conservative and is safe for design. In terms of the influence distance, the 2D analysis agreed with the 3D analysis regarding sections away from the intersecting corner, but the 2D analysis was conservative for the concave slope and was optimistic for the convex slope. Therefore, great care must be taken when interpreting the 2D analysis results when the slope is convex in the longitudinal direction.

Stochastic reliability analysis was conducted to consider both the uncertainty and spatial variability of soils. An example is given herein, where the four material parameters (cohesion, friction angle, particle density, and permeability) are modeled as random fields. Monte-Carlo simulations were conducted to investigate the statistics of the influence and runout distance, associated with the uncertainty of the material properties. Brute-force analysis and machine-learning-aided analysis were compared. For the brute-force analysis, simulations of 1000 samples needed 16 days of computation. However, the estimated probability density function (PDF) of influence and runout distance were still not satisfactory, due to the small sampling size. Therefore, a larger number of samples are required to obtain an accurate PDF. However, traditional stochastic reliability analysis requires a great deal of computing resources and time to achieve a higher accuracy—10⁵ samples require 4–6 years to finish, whereas machine-learning-aided analysis is very efficient.

A neural network was used to establish the relationship between the material properties (which are spatially variable) and the influence and runout distance. With only 500 samples for use as training datasets, the ML algorithms could train a model that predicted the influence and runout distance with good accuracy. With this machine-learning model as a surrogate model to predict the influence and runout distance of millions of samples, a PDF with high accuracy was obtained. Once such a PDF is obtained, a complete stochastic examination can then be conducted, including the mean values of influence and runout distance, identifying the probability that infrastructure lies within the influence of landslides.