Simulating the entire rainfall-induced landslide process using the material point method for unsaturated soil with implicit and explicit formulations

Rainfall-induced landslides cause various types of damage, including damage to infrastructure, leading to devastating economic losses and human casualties. Although various numerical methods have been developed to predict landslide occurrence and the extent of sediment flow, three-dimensional analysis of the entire landslide process in a unified manner is still challenging owing to limitations in computational efficiency and the representation of deformation and flow. In this study, we present a numerical method of rainfall-induced landslides using the coupled hydromechanical material point method (MPM) for unsaturated porous media based on implicit and explicit formulations, which enables the efficient analysis of all rainfall-induced landslide processes, including both quasi-static and dynamic processes. The developed method includes an implicit MPM based on a simplified formulation, which is first applied to the quasi-static analysis in the pre-failure stages in which rainwater infiltrates the ground. The explicit MPM is then applied to the dynamic analysis for post-failure stages in which the ground collapses and flows. A constitutive law for soils is improved in the simulation of landslide initiation and sediment flow by incorporating the effect of cohesion in a visco-plastic model for granular materials. The proposed method was applied to a three-dimensional terrain model of Ashikita town, Kumamoto, Japan, where an actual landslide occurred owing to intense rainfall on July 3–4, 2020. Numerical results from the proposed hydromechanical coupling compared with single-phase MPM indicated that pore water plays an important role in understanding all rainfall-induced landslide processes, from landslide initiation to sediment discharge.


Introduction
Throughout history, landslides have caused severe damage in mountainous areas. Particularly in Japan, landslides occur frequently and cause loss of life and destruction of infrastructure owing to the characteristic terrain and climate of Japan (such as complex terrain with steep mountains, torrential rains caused by rainy season fronts, and typhoons). Owing to the frequent and severe damage caused by landslides, appropriate landslide risk management must be conducted based on more accurate risk evaluations; however, the variation and complexity of the phenomenon lead to difficulty in risk evaluation. Rainfall-induced landslides have a complex mechanism involving the interaction between the ground structure and rainwater, which makes predicting such events particularly challenging.
Among various studies on predicting landslides and reducing landslide damage, numerical methods have been gained increased as a powerful tool to evaluate the complex mechanical behavior of geomaterials and have been applied to landslide disasters. In terms of landslide countermeasures, the following two points must be addressed: landslide initiation points and the extent of the collapsed sediment flow, for which most of the existing numerical methods have been developed.
The limit equilibrium method (LEM) is widely used to evaluate the potential slip surfaces and factor of safety from the balance between shear stress and shear strength for landslide initiation assessments. The LEM can also consider the effect of rainfall when coupled with hydrogeological models that can evaluate the effects of pore water pressure, and various models have been developed in this regard (Montgomery and Dietrich 1994;Pack 1998;Baum et al. 2008). In addition, the effects of the interaction between the soil skeleton and pore water, which often emerge in catastrophic landslides, are not considered by these models. In contrast, the finite element method (FEM) with Lagrangian description was developed as more rigorous method based on continuum mechanics, which enables the detailed evaluation of the stress-strain relation at each point inside the ground. Owing to its representation performance and accurate soil modeling, FEM can predict the processes from the deformation to failure of stable ground. Furthermore, coupled hydromechanical modeling based on the theory of porous media, in which a porous material comprises a solid phase as well as pore fluid phases, has been actively developed for the simulation of various phenomena of saturated/unsaturated soils in FEM (Zienkiewicz et al. 2001;Schrefler and Scotta 2001;Laloui et al. 2003;Ehlers et al. 2004;Uzuoka and Borja 2012). However, these Lagrangian-FEMs are not adequate for post-failure analysis, which exhibits very large deformation and might lead to mesh distortion. Therefore, for sediment flow assessment, computational fluid dynamics (CFD), based on Eulerian description, which observes the solution at fixed points in space, is widely employed (Hungr 1995;McDougall and Hungr 2004;Christen et al. 2010). This framework is suitable for modeling sediment flows as non-Newtonian fluids but not for modeling solids that require history tracking. Thus, it has limitations in terms of pre-failure analysis.
While the above mentioned Lagrangian-and Eulerian-meshbased methods focus only on a partial landslide stage, particle methods have the potential to represent all successive landslides phases because of their meshfree characteristics with Lagrangian descriptions. Numerous attempts have been made to use smooth particle Landslides 20 • (2023) Original Paper hydrodynamics (SPH) (Pastor et al. 2014;Bui and Nguyen 2017), reproducing kernel particle method (RKPM) (Siriaksorn et al. 2018;Wei et al. 2020), material point method (MPM) (Mast et al. 2014;Llano-Serna et al. 2016;Soga et al. 2016;Conte et al. 2019;Xu et al. 2019;Liu et al. 2020;Pan et al. 2021), and others. Among them, MPM, which was originally proposed by Sulsky et al. (1994), has the advantage that its scheme is similar to that of FEM; a continuum body is divided by a set of material points that carry history information, and the field variables are approximated by a background grid. Given that material points and the grid play the equivalent roles of integration points and meshes in FEM, respectively, MPM can be viewed as a variant of Lagrangian FEM in which integration points are independently movable from the mesh. This feature enables users to employ the soil model used in conventional FEM in MPM to represent the post-failure behaviors of structures.
Similar to FEM, coupled hydromechanical MPMs for saturated/ unsaturated soils have also been developed for simulations involving large deformations and flow behavior of soil-water mixtures, such as rainfall-induced landslide and levee failure. Some variants of methods suitable for representing the deformation and flow characteristics of each phenomenon have been proposed, including methods using single-layer material points, which have multi-phase information (Zhang et al. 2009;Higo et al. 2010;Yerro et al. 2015;Bandara et al. 2016), and methods using two-layer material points, each of which represents their own phases (Abe et al. 2014;Bandara and Soga 2015;Yamaguchi et al. 2020).
The above mentioned coupled hydromechanical MPMs can model various types of soil behavior; however, most of their application is limited in the two-dimensional cross-section. Thus, although considering the three-dimensional (3D) effect is vital in simulating rainfall-induced landslides because of its complicated features (such as geometric conditions, rain infiltration, and sediment entrapment), the performance of the method for 3D analysis has been largely unexplored. The limited number of applications of 3D analysis with MPM is caused by its low computational efficiency for simulating all successive processes of rainfall-induced landslides. Specifically, most coupled hydromechanical MPMs are effective for post-failure analysis of landslides because they have been developed primarily as explicit methods for dynamic analysis. However, they are unsuitable for pre-failure analysis, which requires quasi-static analysis of rainwater infiltration over a long period of time, because explicit methods are computationally inefficient for simulating quasi-static processes. Although the majority of MPMs employ explicit approaches, some studies have focused on the accuracy and effectiveness of implicit formulation (Guilkey and Weiss 2003;Sulsky and Kaul 2004). Recently, Zhao and Choo (2020) introduced an implicit scheme in coupled hydromechanical MPM for saturated soil and successfully represented large deformations on a long time scale. Given these studies, it is reasonable to develop implicit schemes for landslide simulation.
In this study, we developed a coupled hydromechanical MPM for unsaturated soils that can represent a series of rainfall-induced landslide processes and demonstrated its representational performance in landslide simulation using a 3D model. To represent both the pre-failure stage, i.e., the quasi-static process of rainfall infiltration, and the post-failure stage, i.e., the dynamic process of sediment flow, we combined the coupled hydromechanical MPMs based on implicit and explicit formulations. Specifically, for pre-failure simulation, we developed an implicit scheme of coupled hydromechanical MPM for unsaturated porous media in reference to coupled hydromechanical FEMs for unsaturated porous media based on a simplified formulation (Ehlers et al. 2004;Uzuoka and Borja 2012), while an explicit scheme based on Bandara et al. (2016) is applied for post-failure simulation. The performance of the proposed method is demonstrated through simulations using a 3D terrain model based on a digital elevation model (DEM) of Ashikita town, Kumamoto, Japan, where a rainfall-induced landslide occurred in July 2020. To demonstrate the effectiveness of the coupled hydromechanical MPM capable of accounting for rainwater infiltration in three dimensions, a single-phase MPM analysis using the same terrain model was also performed for comparison.

Governing equations
Following the theory of porous media for unsaturated mixtures at a finite strain (Schrefler 2002), soils can be described as the superposition of three phases, namely, solid skeleton-water-air, in reference to the solid configuration. For this, we assume the following: • Pore air pressure is negligible (pseudo triphasic formulation).
• Solid grains are incompressible while water is weakly compressible.
• The relative acceleration of pore water to the solid skeleton is ignored.
Unless indicated otherwise, we use subscripts "s" and "w" to represent the solid and water phases, respectively. Then, the momentum balance equations and mass balance equation are written as follows: where = (1 − ) sR + s w wR is the overall density of the mixture, sR , wR is the intrinsic density of solid grain and pore water, is the porosity, s w is the degree of water saturation, u s is the displacement of the solid skeleton, ′ is the Cauchy effective stress tensor, p w is the pore water pressure, 1 is the second-order identity tensor, b is the body force vector, k ws is the permeability coefficient of pore water, g is the acceleration due to gravity, v ws is the relative velocity of pore water, K w is the bulk modules of pore water, and c = − s w p w is the specific water capacity. Herein, we use the over-dot to denote the material time derivatives with respect to the solid skeleton. Equations (1), (2), and (3) describe the momentum balance of the mixture, the momentum balance of the water, and the mass balance of the water, respectively.
Because the relative water acceleration has been omitted, its simplified form can be obtained by substituting Eq. (2) into Eq. (3) as (2) s w wRüs + s w p w − s w wR b + s w wR g k ws s w v ws = 0, (3) s w wR K w − wR c w ̇p w + s w wR ⋅u s + ⋅ s w wR v ws = 0, Landslides 20 • (2023) Equations (1) and (4) are often referred to as simplified formulation or u − p formulation, and this formulation has been used for implicit analysis using FEM in numerous studies (Zienkiewicz et al. 2001;Schrefler and Scotta 2001;Laloui et al. 2003;Ehlers et al. 2004;Uzuoka and Borja 2012) as it requires fewer independent variables. The hydraulic behavior of unsaturated soils is characterized by the soil-water characteristic curve (SWCC), which describes the relationship between the matric suction p c = −p w and the degree of water saturation s w . In this study, the well-known van Genuchten model is employed, which can be written as follows: where s we is the effective degree of saturation, s wr is the minimum degree of saturation, s ws is the maximum degree of saturation, and vg , n vg , and m vg are material parameters.
Whereas Eq. (5) describes the saturation dependency on suction, it ignores the direct effect of porosity on saturation. Although there is a porosity dependency on suction, Eq. (5) sometimes fails to predict suction when the soil skeleton undergoes significant volumetric changes (Aubertin et al. 1998;Gallipoli et al. 2003;Tarantino 2009). Numerical simulations of granular flow, such as those conducted by Dunatunga and Kamrin (2015) and Tampubolon et al. (2017), have indicated that significant dilation may occur in landslides, particularly during the post-failure stage. This highlights the importance of accurately modeling sediment behavior in landslide simulations. Cuomo et al. (2021) also noted that the classical van Genuchten model tends to overestimate suction, which can lead to inaccuracies in predictions of soil displacement and runout distance in simulations using the MPM. They addressed this problem by introducing a maximum limit on suction, whose value is reduced with deviatoric strain accumulation. In this study, to represent the unsaturated state on a broad range of porosities after large deformations, we accounted for the void ratio effect on the van Genuchten model parameter following Tarantino (2009) as follows: where e = ∕(1 − ) is the void ratio and a vg is a material parameter. Equation (6) allows the representation of a decrease/increase in suction due to an increase/decrease in porosity and prevents suction over-/under-estimation in landslide analysis.
In the van Genuchten model, the water permeability is given as follows: where k s ws is the maximum value of water permeability, and vg and m vg are material parameters.
The solid skeleton is modeled as a visco-plastic model for granular deformation based on the finite strain theory. Employing the numerical algorithm established in Simo and Taylor (1985) and Simo (1992), (5) s w = s wr + (s ws − s wr ) 1 + ( vg p c ) n vg −m vg , , s we = s w − s wr s ws − s wr , we assume the following multiplicative decomposition of the deformation gradient tensor: Here, F e and F p are the elastic and plastic components, respectively. To represent the elastic deformation, we employ the Hencky's hyperelastic law (Simo 1998). Although more sophisticated elasticity models for granular materials could improve the representation performance, the simplicity of Hencky's model is advantageous in terms of computational efficiency.
Recently, the visco-plastic model proposed by Kamrin (2010), which is based on a rheology model for dense granular materials (Jop et al. 2006), was applied in some MPM studies for representing the deformation and flow behavior of granular media (Dunatunga and Kamrin 2015;Xu et al. 2019). In this study, we followed these ideas and added the cohesion effect to represent landslides that have a transition process from a soil structure to dense granular flow with softening. The elastic limit is represented by the Drucker-Prager like yield criterion as follows: with where is the Kirchhoff stress tensor, J 2 (s( )) is the second invariant of the deviatoric stress tensor s = − 1 3 tr( )I , p = 1 3 tr( ) is the hydrostatic pressure component, c is the cohesion, and s is the static internal friction angle. Because the associated flow rule often generates excess inelastic dilation, a non-associated flow rule is defined as follows: with where I 0 is a dimensionless constant, d is the soil particle diameter, and 2 is a limiting values of the internal friction angle. With the piecewise linear hardening/softening assumption, the evolution of cohesion with accumulated plastic deformation ̄p can be described as follows: where H s is the hardening/softening parameter. Here, the lower limit of cohesion is set considering the softening behavior of soil in the large deformation. Here, as the plastic potential Ψ contains only a deviatoric component, plastic deformation is assumed to be incompressible when the plastic yielding occurs at the cone surface of the yield surface. However, if the trial stress projection is outside the cone surface, the stress state is projected onto the apex of the Landslides 20 • (2023) Original Paper yield surface (de Souza Neto et al. 2011). This leads to the inevitability of excess volumetric plastic dilation under the framework of classical elasto-plastic theory. To overcome such issues, some studies (Dunatunga andKamrin 2015, 2017;Hurley and Andrade 2017;Tampubolon et al. 2017) have corrected the volume component in elastic deformation. Based on these ideas, to evaluate the trial stress state, we define the modified elastic deformation gradient as follows: where the superscript " tri " denotes the variables at the trial state. Equation (14) indicates that the volumetric component of the trial elastic deformation is evaluated by the total deformation and is not affected by the previous state of the plastic deformation.

Implicit formulation of the coupled hydromechanical MPM
The pre-failure stage of the landslide is simulated by an implicit MPM for quasi-static analysis based on the u − p formulation. The strong form of the u − p formulation consists of Eqs. (1) and (4) and the following boundary conditions: where ū s is the prescribed displacement of the solid skeleton, t is the prescribed traction, q ws is the prescribed flux of pore water, p w is the prescribed water pressure, and n is the outward unit normal vector on the boundary. The boundary surfaces hold Ω = Ω s,u ∪ Ω s,t = Ω s,q ∪ Ω s,p and � = Ω s,u ∩ Ω s,t = Ω s,q ∩ Ω s,p . Then, the weak form of the u − p formulation for quasi-static analysis can be obtained via the following standard procedures: where v s and p w are test functions. w s and w p are the weak form of the momentum balance of the mixture and the mass balance of the water, respectively. In addition to Eqs. (19) and (20), boundary terms are added to realize the Dirichlet boundary conditions (see Appendix 1). The MPM discretization is performed by the set of Lagrangian particles, which represents the solid skeleton body and the computational grid nodes where the discretized equations are solved.
We use subscripts "p" and "I" to denote the variables related to the particles and grid nodes, respectively. We assume that the physical domain is divided into a set of n mp subdomains, and the volume of each subdomain is concentrated at the corresponding single material point, implying that the integrals for arbitrary field appearing in the weak form can be evaluated as follows: where x mp and V mp are the position vector and the point volume associated with a material point mp , respectively. In a manner similar to the physical domain, the boundary domain is also divided into a set of n bp boundary points and the boundary integrals are evaluated as follows: where A bp is the point area associated with a boundary point bp.
For spatial discretization, the continuum fields are approximated as follows: where n I is the number of grid nodes, B I (x) is the nodal basis function evaluated at x , and I is its coefficient. There are several choices for the basis function for MPM; we employ the extended B-spline basis for implicit formulation (see Appendix 2). For time discretization, the central difference scheme is employed for both the displacement and pressure field. Discretization of Eqs. (19) and (20) in space and time presents discrete residual vectors n+1 s and n+1 p with nodal components defined as follows: where B Ip is the basis function for node I evaluated at the position of each particle p, A bp = A bp n bp is the area vector, and x n+1 denotes the spatial gradient at time step n + 1 . In addition to the above terms, stabilization terms are employed for pressure stabilization (see Appendix 3).
For the non-linear system solution, we apply the iterative procedure along with the Newton-Raphson method, and the following linearized equations are solved in every iterative step k + 1: Here, d s and d w are the corrector vectors of solid displacement and pore water pressure, respectively. Details regarding the tangent matrices can be found in Uzuoka and Borja (2012).

Explicit formulation of the coupled hydromechanical MPM
For the dynamic analysis of the post-failure stage of landslides, an explicit MPM based on the u − w formulation is employed in reference to Bandara et al. (2016). Similar to Eqs. (19) and (20), the weak form of the u − w formulation is obtained from the linear momentum balance (Eqs. (1) and (2)) and boundary conditions (Eqs. (16) and (18)) as where v ws is the test function. w s and w w are the weak form of the momentum balance of the mixture and water, respectively. The Dirichlet boundary condition is then obtained by adding a boundary term, as in the implicit solution method (see Appendix 1). For time discretization, the forward Euler method is applied for solid displacement fields. Application of Eqs. (21)-(22) and the approximation of velocities and test functions of the solid and water phases in Eqs. (1) and (2) with Eq. (23) yield the following discretized equations: Here, the nodal components of each term are defined as n̈ n where IJ is a Kronecker delta function.

Numerical procedure
The numerical simulation of rainfall-induced landslides can be divided into two stages: the quasi-static and dynamic stages; in other words, these are the stages before and after slope failure, respectively. In the quasi-static stage, implicit MPM is applied with prescribed rainfall boundary conditions, in which the time increment is controlled adaptively in reference to Ehlers et al. (2004). The numerical algorithm of an implicit MPM is written as follows: 1. Restore the grid. 2. Compute tangent stiffness matrix and residual vector in Eq. (26). 3. Solve linearized equation (26) for iteration k + 1 and update solid displacement increments and pore water pressure increments on the grid: 4. Update state variables at material/boundary points: Landslides 20 • (2023)

Original Paper
5. Update Cauchy stress using the constitutive functions. 6. Check convergence for equilibrium: If convergence has not been achieved, return to 2. 7. Save state variables at material/boundary points. 8. Update spatial coordinates of material/boundary points: If convergence of the Newton-Raphson method is not attained within the specified time increment or the prescribed rainfall time is elapsed, the method is switched to explicit MPM, and the states of the particles in the last convergent state of the implicit MPM are passed on. The numerical algorithm of an explicit MPM is written as follows: 1. Restore the grid and map the velocity of material points to grid nodes: 2. Compute the solid acceleration vector ü n s, I by Eq. (29) and update solid velocity on the grid: 3. Compute average relative water velocity w n ws, I by Eq. (30). 4. Update state variables at material/boundary points by Eqs.
(50) u n+1 s, p =u n s, p + Δt  x n B Ip ⋅ w n ws, I .
In addition to the above basic algorithm, we utilize the MUSL procedure (Sulsky et al. 1995) and pressure stabilization method (see Appendix 3) for robustness improvement. The numerical experiment reveals that the method can reproduce the fundamental behavior of unsaturated porous media (Appendix 4).

Application to an actual landslide event
Simulation setting A landslide triggered by intense rainfall on July 3-4, 2020, in Tagawa, Ashikita town, Kumamoto Prefecture, Japan (Fig. 1), was selected to demonstrate the effect of the 3D coupled hydromechanical MPM. The test site is approximately 400 × 300 m and is composed of mountainous terrain with an elevation ranging from 10 to 150 m. Under normal conditions, the collapsed slope was a shallow catchment area without flowing water. According to Jitousono et al. (2021), the maximum and average collapse depths were approximately 10 and 3 m, respectively, which are more significant than those in a typical landslide. Figure 2 shows the time series of hourly rainfall observed in the Ashikita observation station (Jitousono et al. 2021). The total amount of precipitation between the beginning of the rainfall event (4:00, July 3) to landslide occurrence (5:00, July 4) was 384 mm. The hourly rainfall from 4:00 to 5:00 when the landslide occurred was the highest during the period, exceeding 80 mm.
For the initial setting of the ground geometry for numerical analysis, 5 m mesh DEM data published by Geospatial Information Authority of Japan (GSI) website (2022) were used. Figure 3 shows the topological information of the target area on a plane view, including (52) x n + 1 p = x n p + Δt In this study, the ground model is assumed to comprise a single soil layer overlying rigid, impermeable bedrock. Thus, the depth of the topsoil (i.e., the relative position of the bottom surface to the top surface of the soil layer) corresponding to the DEM is set at first. The soil depth is known as an important factor affecting the stability of a slope and is usually measured via penetration tests. Some penetration test results are available in KuniJiban (2022) in Japan; however, these points are unevenly distributed and inadequate to evaluate the landslide risk over a wide area. Some previous studies (Sarkar et al. 2013;Salciarini et al. 2017) employed the kriging method to resolve this issue, which predicts the continuous data distribution from discrete spatial data. Another approach based on physical terrain features can be found; a relationship between soil depth and the ground surface profile, such as slope angle and curvature, can be evaluated over multiple regions (Heimsath et al. 1997;Tesfa et al. 2009;Patton et al. 2018). Given the lack of data in numerous regions, the latter modeling method was employed in this study to investigate the influence of the parameters. Patton et al. (2018) proposed the following equation to generalize the relationship between soil depth H and slope curvature C: where ΔH ΔC is the slope of H -C function and H is the average of the depth within a target region. Some methods have also been proposed to estimate ΔH ΔC and H ; however, their applicability to arbitrary regions, especially to characteristic terrains such as the mountainous regions of Japan, has not been well verified. Therefore, we tested two settings based on a field survey of target landslides to investigate the effect of soil depth distribution: one was the curvature-independent soil depth with H = 6 m, i.e., uniform soil depth, and the other was the curvature-dependent soil depth with H = 50C + 6 m.
Subsequently, the surfaces of the top and bottom surfaces of the soil layer can be generated by the linear interpolation of the DEM and the corresponding topsoil depth, respectively. Prior to the generation of the particle model along the terrain, particles were tentatively placed in the entire computation domain. We utilized the uniform Cartesian grid for computation with cell length h and placed six material and boundary points per h. Next, the material points located in the topsoil layer, i.e., the layer between the top and bottom surfaces, were selected. Furthermore, the boundary points located in the nearest neighborhood of the boundary surfaces were extracted. Figure 4 shows the numerical model with the curvature-dependent soil depth.
The material properties of the soil are provided in reference to the penetration test data near the site obtained in KuniJiban (2022) and the typical values based on the soil classification. Note that although these parameters should ideally be identified based on in-situ measurements and parameter studies, as the purpose of this analysis is to demonstrate the performance of the coupled hydromechanical MPM rather than to accurately reproduce actual conditions, the parameters are determined within a range that can be estimated from publicly available information. According to penetration test results, the type of topsoil in this area at a depth of approximately 6 m below the ground surface was classified as sandy gravel based on the soil classification in Japan, which is roughly equivalent to well-graded gravel according to the US Geological Survey (USGS). Based on the material property values for each soil type investigated in previous studies (e.g., porosity (Das and Das 2008), friction angles (Koloski et al. 1989; Geotechdata Info 2013), empirical constant I 0 (Jop et al. 2006;Lucas et al. 2014), van Genuchten parameters (Wösten et al. 1999), and permeability (Carter and Bentley 1991)), the material parameters listed in Table 1 are used on the entire body of the soil layer models in all simulation cases presented in this section. Here, the cohesion is set taking into account that the N value in the penetration test is approximately 5 to 10 and that vegetation is present (Ali and Osman 2008).
In all stages, the cross-sections of the soil layers were fixed on the horizontal direction and drainage, while the bottom surface of the soil layer was rigidly fixed and impermeable. Furthermore, the top surface of the soil layer was free to move, and the rainfall boundary condition was uniformly applied as the prescribed mass flux obtained via linear interpolation of the rainfall history from 4:00 a.m. on July 3 to 6:00 a.m. on July 4 in Fig. 2 in the pre-failure simulation with implicit MPM. Prior to applying the above condition, a steady state was created by applying a constant rainfall condition of 5 mm/h for three days, based on the normal monthly rainfall in the Kumamoto Prefecture during June and July, which is approximately 400 mm. In the following explicit MPM post-failure simulation, because the effect of rainfall during the tens of seconds on sediment flow was small, the rainfall boundary condition was omitted, and the pore water flux at the topsoil surface was assumed to be zero.
Because in MPM, the computational grid is independent of the particles and is reset at each step, the grids can be easily reconstructed during computation. As the computation time for the implicit method is largely dependent on the number of degrees of freedom of the computational grid compared with the explicit method, and a finer grid is desirable for sediment flow analysis,

Original Paper
the grid was changed in quasi-static and dynamic analysis. Specifically, grids of linear extended B-splines with h = 2m and quadratic B-splines with h = 1m were used for the quasi-static and dynamic analyses, respectively. The time increment Δt in the implicit MPM was initially set to 1200 s but was further subdivided with a lower limit of 1 s if convergence could not be achieved. In explicit MPM, the time increment was constant and set to 10 −4 s.
In addition to the coupled hydromechanical analysis, to investigate the effect of pore water, we conducted single-phase analyses for the soil skeleton only in which the terms relevant to the water phase were omitted in the aforementioned implicit and explicit MPMs. In this case, the landslide was triggered by the strength reduction method (Matsui and San 1992;Dawson et al. 1999), i.e., the cohesion of the entire body was uniformly and gradually reduced from its initial value to zero until the Newton-Raphson method no longer converged in the quasi-static analysis. As in the coupled analysis, the cohesion reduction value in the Newton-Raphson was adjusted between 1/3 and 1/3072 kPa, and the value at which it no longer converged was used for the following dynamic analysis. All other parameters and numerical setting were the same as those in the coupled hydromechanical analysis. Figure 5 shows the evolution of pore water pressure on the bottom surface of the soil layer obtained by quasi-static analysis using two numerical models with different soil depths. Furthermore, the evolution of the accumulated plastic strain is shown in Fig. 6, where only the parts with ̄> 0.01 are shown. In both cases, solutions were obtained without lack of convergence in the analysis under rainfall conditions up to 6:00 a.m. on July 4. The pore water pressure distribution after 4:00 a.m. on July 4 shows positive values over a wide area of the target region, along with a corresponding increase in the accumulated plastic strain. However, there were differences in their distribution; i.e., the plastic strain shows large values mainly along the areas marked as A to D in Fig. 6, which have a steep slope, while the pore water pressure is concentrated in valleys. A comparison of the results of the two models indicates

Original Paper
that the distribution of each value is generally similar, but the case with the curvature-dependent soil depth shows relatively higher pore water pressure values in the valley. This is a reasonable result because the slope of the valley topography on the bottom surface is larger in the model with curvature-dependent soil depth. In the final state after prescribing rainfall, the bottom surface of the soil layer shows a positive pore pressure value, i.e., it is saturated, while the top surface of the soil layer is partially saturated, as shown in Fig. 7. Although a slight local suction distribution due to progress of the solid deformation can be found, the saturated areas are mainly in valleys, plains, and the toes of slopes, and we can see that the case with the curvature-dependent soil depth has a larger saturated area around valleys than the case with the uniform soil depth. The results of the single-phase analysis with different soil depths are presented in Fig. 8, where only the material points having ̄> 0.01 are picked in the state just before the slope becomes unstable. Similar to the coupled analysis, plastic deformation occurs mainly on steep slopes, but is moderate in the vicinity of area A, where the landslide actually occurred.
Given that the material parameters of the solid structure in the single-phase analysis are the same as in the coupled analysis, these differences in deformation distribution indicate the effect of rainwater infiltration. As water pressure contributes to the total stress, as seen in Eq. (1), suction works essentially in the same way as cohesion, and the reduction in suction or cohesion leads to a reduction in the shear strength of the soil skeleton due to its pressure-dependent visco-plastic model. In the single-phase analysis, plastic deformation is caused by a uniform reduction in the overall cohesion. In contrast, in the coupled analysis, plastic deformation is caused by a non-uniform reduction in suction due to the infiltrated rainwater 3D distribution. As a result, in the coupled analysis, the plastic deformation seems to be relatively larger in the valley topography, where the infiltrated rainwater is concentrated, compared with the single-phase analysis. These results suggest that a coupled analysis capable of incorporating the 3D seepage water behavior might be able to capture the failure modes of landslides that cannot be captured by a single-phase analysis with the strength reduction technique.
The dynamic analysis results with coupled hydromechanical MPM and single-phase MPM are shown in Figs. 9 and 10, respectively. Here, material points with accumulated plastic strain ̄> 1 are shown. Each of these calculations with the coupled hydromechanical MPM and single-phase MPM took approximately 48 and 24 h, respectively, using 480 CPU cores.
We can see that the onset points of sediment flow in the coupled hydromechanical MPM results do not necessarily coincide with the locations where relatively higher plastic deformation emerges in the quasi-static analysis (indicated by A to D in Fig. 6). Specifically, sediment flows occur only around areas A and C in the result with the uniform soil depth, and around areas A, B, and C in the result with the curvature-dependent soil depth. The difference in the sediment flow around area C between the two models with varying soil depths may be related to the deeper soil layer and higher pore water pressure in the model with the curvature-dependent soil depth because area C has a valley topography. Meanwhile, no sediment flow occurred in any of the results in area D, which has a steeper slope but a lower pore water pressure than the other areas in the figure.
In contrast, the results of the single-phase analysis exhibit a different trend, i.e., the regions with relatively large plastic strain in the quasi-static analysis (Fig. 8)

Original Paper
the landslide actually occurred. When compared with the results of coupled hydromechanical analysis, a significant difference can be observed in the flow distance, with results of single-layer analysis exhibiting a shorter flow distance. Figures 11 and 12 show the distribution of debris flow on a planer map and the surface profile of the soil layer at 30 s on a crosssection of a valley where a landslide actually occurred. The surface profiles of the soil layers are represented by sets of points at 5 m intervals, where each point elevation value is the coordinate of the material point located at the highest position within each interval. The coupled hydromechanical MPM results show that there is slope failure and erosion as well as sedimentation on the plains, which is close to the actual situation investigated by Jitousono et al. (2021). In contrast, the single-layer MPM results show sedimentation in the mid-slope due to sediment inflow from the adjacent ridge but no slope failure.
The difference between the single-phase MPM and coupled hydromechanical MPM in the dynamic analysis suggests that, in addition to the difference in plastic deformation distribution in the pre-failure analysis, the influence of pore water is also significant during large deformation and flow. In the single-layer MPM, the cohesion decreases with the development of plastic deformation due to the plastic-softening constitutive law. Conversely, in the coupled hydromechanical MPM, the pore water pressure changes as the soil skeleton volume changes in addition to plastic softening. As a result, the change in water pressure due to the partial volume change causes a change in shear strength, possibly resulting in a significant difference in post-failure behavior.
These numerical results demonstrated the possibility of the coupled hydromechanical MPM considering the effect of pore water to represent the mechanism of the rainfall-induced landslide, which cannot be achieved with single-phase MPM. In particular, the coupled hydromechanical analysis is required to represent landslides on valley topography where landslides actually occurred in the model used in this study, which assumes a uniform distribution of soil types. However, landslides were observed in areas where they did not actually occur in both numerical models with different soil depth distributions. Additionally, the resulting collapsed soil accumulation caused longer distances of sediment flow than the actual situation. These discrepancies can arise from the lack of consideration of actual non-uniform soil property distribution, and detailed modeling of the ground may reproduce the phenomena accurately. The parametric study demonstrates that different soil properties can lead to different failure modes (Appendix 5). Nonetheless, the coupled hydromechanical MPM still provides useful insights into landslide risk evaluation that cannot be obtained from conventional methods.

Conclusion
In this study, we developed a numerical method for simulating all stages of rainfall-induced landslides and demonstrated the representation performance of the proposed method through 3D simulations of an actual rainfall-induced landslide event in Japan. The proposed method results were compared with single-phase analysis results. A coupled hydromechanical implicit MPM for unsaturated soil was constructed based on a simplified formulation and was used for rainfall infiltration analysis in the pre-failure stage of landslides, while a coupled hydromechanical explicit MPM for unsaturated soil was applied for sediment flow simulation in the post-failure stage. The combination of implicit and explicit MPMs enabled us to simulate a series of landslides, including quasi-static and dynamic processes, using consistently visco-plastic models for granular materials with cohesion.
The proposed method was applied on a 3D terrain model of Tagawa, Ashikita town, Japan, using rainfall boundary conditions based on observational data. A comparison between the proposed coupled hydromechanical and single-phase MPMs showed that pore water significantly affects the 3D analysis, resulting in changes in the location of landslide initiation and the flow distance of collapsed sediment. Furthermore, we found that variations in the distribution of soil layer depth affected the pore water distribution and influenced the landslides occurrence. These results prove that the developed coupled hydromechanical MPM with a 3D model is an effective numerical method for reproducing and investigating the complex landslide mechanism in detail. It is worth noting that the constitutive model used in this study is built upon previous efforts, however, there is potential for further improvement in order to more accurately model the mechanisms of rainfall-induced landslides. Further investigations using the developed method with more detailed data (e.g., soil material parameters and rainfall intensity) and applications in different

Original Paper
where V c is the grid cell volume. The stabilization parameter is set to = 10 3 wR ∕ 2GΔt in this study. Furthermore, although not as serious as the liquid phase, the solid phase can also produce inappropriate results owing to its pressure-sensitive plasticity model described in "Governing equations" section. Thus, we employed the F-bar method for implicit MPM (Coombs et al. 2018), and the deformation gradient was modified as follows: where ΔF 0 is the deformation gradient increment evaluated at the geometric center of material points within a cell x 0 where the material point mp is located.
In explicit MPM, pressure oscillations of pore water caused by a high volumetric modulus are known to lead to inaccurate numerical results and should be addressed along with the solid phase as in implicit MPM. Thus, we combined the smoothing method (Mast et al. 2012) and reduced integration method (Bandara and Soga 2015) for the pore water pressure that suffered from severe instability. Smoothing of the arbitrary field evaluated at the material points is defined as follows: Equation (60) is applied to the increments of pore water pressure Δp w and the determinant of the deformation gradient det(ΔF) . Then, the modified version of the deformation gradient is obtained as follows:

Appendix 4 Leaking flow analysis
We reproduce the experiment of leaking flow in a soil column, as conducted in Liakopoulos (1964), to evaluate the performance of the proposed coupled hydromechanical material point method (MPM). This problem has been widely used as a benchmark in various studies (Schrefler and Scotta 2001;Laloui et al. 2003;Ehlers et al. 2004;Uzuoka and Borja 2012;Bandara et al. 2016) for modeling unsaturated porous media. In the experiment, soil columns were filled into impermeable cylindrical vessels with open bottoms and tops, allowing for one-dimensional flow and of the water outflow velocity at the bottom deformation to occur.
Our numerical model simulates a 1 m tall soil column divided into 20 uniform grid cells, with each cell containing 4 material points in the vertical direction. We use linear B-splines and impose boundary conditions at grid nodes. The bottom surface is fixed in the vertical direction while the top surface is allowed to move freely. Water in/out is initially allowed on both the top and bottom surfaces, but later only on the bottom surface in the leaking Minimum degree of water saturation s wr 0.02 Maximum degree of water saturation s ws 1.0 (a) (b) Fig. 13 Numerical results of leaking flow analysis: a distribution of the pore water pressure and b time history of water outflow velocity process. The horizontal movements of all phases are fixed on all nodes. Assuming the solid phase to be Hencky's hyperelastic material, material parameters are set as shown in Table 2, based on previous numerical studies (Ehlers et al. 2004;Uzuoka and Borja 2012;Bandara et al. 2016). The simulation begins with a stationary state obtained through a quasi-static analysis using the implicit MPM. The leaking process is simulated for the first 3 minutes using the implicit MPM with a time step of Δt = 10 −1 s. After this, two cases are compared: the case using the implicit MPM and the case switching to the explicit MPM with a time step of Δt = 10 −4 s.

Original Paper
Figure 13(a) shows the development of the pore water pressure distributions in the implicit and explicit MPMs. Additionally, Fig. 13(b) shows the time history of the water outflow velocity at the bottom material point. Although the slower outflow tendency, which is consistent with the previous FE analysis (Ehlers et al. 2004), can be observed, the results obtained using the implicit and explicit MPMs are in good agreement with each other, and they accurately reproduce the trend of the time series of outflow in the experimental results. This suggests that the proposed method has the capability to capture the fundamental mechanics of unsaturated porous media.

Effect of material parameters
To understand how material parameters affect landslide failure modes, we conducted six simulations with different parameter values: minimum friction angle ( s ) of 38 • and 28 • , initial cohesion ( c 0 ) of 30 kPa and 10 kPa, and saturated permeability coefficient ( k 0 ws ) of 10 −3 m/s and 2.5 × 10 −4 m/s. The default case is a uniform soil depth model with the same initial and boundary conditions outlined in the "Application to an actual landslide event'' section. In these simulations, all other parameter values are fixed as in the default case. Figure 14 illustrates the distributions of accumulated plastic strain ( ̄p ) at the final state in the quasi-static analysis of each case, where only material points with ̄p > 0.01 are displayed. It is noteworthy that the case with c 0 = 10 kPa did not yield a convergent solution after 10 p.m. on July 3, while the other cases yielded convergent solutions in the analysis under rainfall conditions up to 6:00 a.m. on July 4. Figures 15 and 16 show the distributions of pore water pressure on the bottom surface of the soil layer and the suction distributions on the top surface of the soil layer in the cases with c 0 = 10 kPa, k 0 ws = 10 −3 m/s, and k 0 ws = 2.5 × 10 −4 m/s, respectively. The cases with c 0 = 30 kPa, s = 28 • , and s = 38 • are not shown as they show similar results to the default case in Figs. 5(a) and 7(a). The results indicate that the trend of plastic deformation distribution is consistent across all cases, with deformation concentrated mainly in areas A to D in Fig. 14, similar to the default case in Fig. 6(a). However, the amount of deformation varies. The cases with s = 38 • and c 0 = 30 kPa, which have high material strength, exhibit less plastic deformation than the case with s = 28 • , which has lower material strength. Even though the case with c 0 = 10 kPa received less rainfall than the other cases, resulting in lower water pressure, it still displayed large plastic deformations mainly in area D. The comparison of the water pressure between the cases with k 0 ws = 2.5 × 10 −4 m/s and k 0 ws = 10 −3 m/s in Fig. 15 shows that the former has a higher water pressure over a larger area, specifically in areas with valley and plain topography. Additionally, the saturation distribution in Fig. 16 indicates that the case with k 0 ws = 2.5 × 10 −4 m/s has a wider area of saturation. As a result, it can be consistently observed that the case with smaller permeability leads to greater plastic deformation. Figures 17 and 18 present the results of the dynamic analysis in the post-failure stage, where only material points with the accumulated plastic strain of ̄p > 1 are displayed. It is worth noting that the cases with c 0 = 30 kPa and k 0 ws = 10 −3 m/s have been excluded from the figures as no landslides occurred in these cases. The figures clearly demonstrate the differences in the onset points and flow paths of the sediment flow for each case. In the case with s = 38 • , sediment flow is observed only in the vicinity of area C. On the other hand, in the case with s = 28 • , landslides occur not only in areas A and C, as in the default case in Fig. 9, but also in area B. The variation in the landslide occurrence areas in the above cases can be attributed to the difference in material strength. Furthermore, in the case with c 0 = 10 kPa, landslides occur exclusively in area D, which is not the case in the default scenario. The findings of the study indicate that the location and behavior of landslides can vary based on the amount and distribution of pore water in the soil, depending on the rainfall history. The simulation results show that the case with k 0 ws = 2.5 × 10 −4 m/s resulted in landslides in four areas (A to D) that were more extensive than the default case. In contrast, the case with k 0 ws = 10 −3 m/s did not result in any landslides, suggesting that landslides are more likely to occur with smaller water permeability. Additionally, the cases with s = 28 •

Original Paper
and k 0 ws = 2.5 × 10 −4 m/s, where the landslide area was larger than in the default case, resulted in faster sediment flow and earlier reach of the plain area. Overall, the results suggest that higher material strength and lower water permeability result in larger landslide and sediment flow areas or earlier landslide occurrence.