Sensitivity and Calibration of Three-Dimensional SPH Formulations in Large-Scale Landslide Modeling

Numerical prediction of landslide runout and deposition is important for estimating landslide risk and developing mitigation plans. The choice of a suitable model and its parameters and a confident calibration strategy are crucial for numerical simulations. Here, we evaluated two constitutive models with a three-dimensional smoothed particle hydrodynamics (SPH) method by simulating the catastrophic 11 October 2018 Baige landslide. The results indicate that both the soil mechanic and fluid models can capture the dynamic runout and deposition morphology while using different values of input parameters. A point-wise comparison of deposit elevation can minimize the calibration error. Numerical models were constrained accurately by utilizing both the static observation data and dynamic seismic signals. The effects of friction on deep-seated landslides motion

Landslides often occur in complicated terrain where the sliding material includes forward motion and lateral spreading in response to local topography.Hence, developing three-dimensional (3D) models is crucial to describe the dynamical and internal stress of landslide motion in real-world terrain (Aaron, 2017;Dai et al., 2014).However, the 3D models for large-scale simulations usually require a large amount of computation time.One solution to address this issue is introducing the graphics processing unit (GPU) acceleration method (Peng et al., 2019(Peng et al., , 2022)), which is a type of robust, cost-and energy-effective computing device that has been popular in numerical modeling in recent years.Currently, a single high-end GPU has thousands of processing cores that are capable of simulating tens of millions of SPH particles in parallel (Crespo et al., 2015;Peng et al., 2019Peng et al., , 2022;;Zhan et al., 2019), thus, enabling researchers to conduct 3D SPH simulations of natural landslides.
The first practical challenge in 3D SPH simulations is to establish an appropriate constitutive relationship to describe the behavior of the material (McDougall, 2017).Generally, as a first approximation, landslide mass is usually modeled as single-phase homogeneous material either using soil mechanic models or non-Newtonian fluid models (Aaron, 2017;Huang & Dai, 2014;Hungr & McDougall, 2009;Pastor et al., 2018).However, whether the soil or fluid model is more suitable for simulating high-speed landslides has not been thoroughly investigated.Second, the mixtures in landslides mass are usually assumed as Mohr-Coulomb materials, thus, the shear strength (friction and cohesion) significantly influences the movement of landslides (Yamada et al., 2018) and their scale (Jeandet et al., 2019).However, there is a lack of sensitivity analysis of these two parameters in constraining landslide movements.Although some researchers (Aaron et al., 2018(Aaron et al., , 2019;;Yamada et al., 2013Yamada et al., , 2018) ) have explored this issue, they mainly conducted in the depth-averaged models.The velocity in the depth-wise direction is neglected and the pressure of the material is considered to be hydrostatic in the assumptions of the depth-averaged models (e.g., Gao et al., 2017;Hungr & McDougall, 2009;McDougall & Hungr, 2004;Ouyang et al., 2019;Pastor et al., 2009Pastor et al., , 2014;;Savage & Hutter, 1991).This may not work very well on steep slopes as many depth-averaged models also use a bed-normal approximation, which means on steep slopes there will be a significant vertical velocity considering Cartesian coordinates.In principle, real 3D models can overcome these limitations, yet, it has not been performed so far.Third, the calibration of numerical simulation still needs to be improved.Previous researches, for instance, Dai et al. (2014), Ouyang et al. (2019), and An et al. (2021) mainly focused on comparing numerically and field-obtained deposition data (e.g., area, runout distance, and maximum height), which can be conveniently derived from post-event field investigations.In some cases, this simple calibration can help us to estimate the bulk characteristic of landslides.Other researchers used the forced vortex equation (Aaron et al., 2019) and dynamical method to estimate superelevation velocities (Pudasaini & Jaboyedoff, 2020), or building damage to estimate flow velocities (Kean et al., 2019).However, the post-event field measurements cannot be used directly to provide velocity-time and force-time estimates of the sliding process.
Alternatively, seismic data provides a unique opportunity to study the time-dependent landslide dynamic processes (e.g., Brodsky et al., 2003;Ekström & Stark, 2013;Hibert et al., 2014).Seismic signals are generated by the interaction between the flowing mass and the bedrock during the initialization, acceleration to deceleration of the sliding mass (W.Li, Chen, et al., 2019).The analysis of landslide-generated seismic waves can complement those from post-event field investigations.The long-period seismic waves (a few seconds to hundreds of seconds) generated by landslides can be recorded by seismometers over distances of tens to hundreds of kilometers and then be inverted to estimate the interaction force between the sliding body and the Earth's surface over time (Moretti et al., 2012).The forces and resulting seismic signals contain information about the topography, the mass flow dynamics (e.g., velocity, volume, and density), and the mechanical properties of the landslide deposits (e.g., friction coefficient) (Moretti et al., 2015;Yamada et al., 2018).By combining seismic wave analysis and post-event field observations, we can solidly calibrate models and infer the landslide processes based on the mass, momentum, and trajectory of landslides (Allstadt, 2013;An et al., 2021;Favreau et al., 2010;W. Li, Chen, et al., 2019;Moretti et al., 2012Moretti et al., , 2015;;Zhang et al., 2019).
To give a practical guide on numerical modeling for geophysical flows, we investigated two rheological relations with SPH formulation applied to the 11 October 2018 Baige landslide event and explored the inclusion of both seismic data and field observations into the calibration process.Moreover, some indicators were used to check the goodness of the simulations.In the remaining part of the paper, we first present the background of the Baige landslide in Section 2. In Section 3, the numerical method and the analysis approach for obtaining time-force curves by inverting the seismic data are briefly introduced.Next, we present the sensitivity analysis of two constitutive models and the dynamic processes of the Baige landslide (Section 4).In Section 5, we discuss the performance of the models, model parameters, shear strength properties inferred from the Baige landslide, and model limitations.

Study Site and the Baige Landslide
On October 11th and 3 November 2018, two successive landslides occurred at the Baige village (31° 4′ 56.41′′ N, 98° 42′ 17.98′′), on the border between the Jiangda county and the Baiyu county, China (see Figure 1a).In this study, we focus on the first event.The landslide was located in a tectonically active and highly fractured geological region in the southeastern part of the Qinghai-Tibet Plateau.The region is mainly composed of Triassic and Proterozoic strata.The geological setting is controlled by folds and faults in the north-western direction.A total of 6-7 trunk faults, most of which are median-angle thrusts with right-lateral slip, developed in the landslide area.One of the Boluo-Muxie faults crosses through the back part of the Baige landslide (Cai et al., 2020).The upper part of the landslide, namely the fault hanging wall, is the T 2−3d -dk stratum, and the lithology is gneiss and greenschist.The lower part of the landslide, that is, the footwall of the fault, is the J ηγ stratum, and the lithology is granite (Cui et al., 2020).A 300-m-wide band of severely fragmented rock (formed by regional faulting and deep weathering) composed of plagioclase gneiss and serpentinite, which makes it prone to sliding (Fan et al., 2019).The dip direction of the bedding planes and foliation is about 230°-250° with a dip angle of 40°, which is reversed with the dip direction of the slope (75°) (Cai et al., 2020).
The landslide site has a typical continental monsoon climate, where the annual rainfall is unevenly distributed, such that almost 78% of precipitation occurs during the summer (H.-B.Li, Qi, et al., 2019).According to the nearby Dege Country Meteorological Station (about 20 km away from the site), the region has no significant rainfall event before 11 October 2018 (Tian et al., 2020).In October 2018, the maximum daily rainfall was 16.1 mm, typical for the rainy season at the site.Five daily rainfall events exceeded 16.1 mm in August and September 2018 (see Figure 2a).The accumulated annual precipitation in 2018 before the landslide event reached 790.4 mm, the highest total precipitation during the past seven years and higher than the annual rainfall of 2016 (606.6 mm) and 2017 (702.7 mm, see Figure 2b).This indicates that the Baige landslide may have been triggered by accumulated rain over a long period of time rather than a single heavy rainfall.The long-term rainfall infiltration might have caused a progressive and persisting build-up of hydrostatic pressures within the fractures, possibly promoting slope slow creep deformations at the boundaries and within the source mass (Fan et al., 2017).This led to crack opening and propagation over a long time, which finally accelerated toward a dramatic mass collapse (Fan et al., 2019;H.-B. Li, Li, et al., 2019).This speculation was supported by multitemporal remote sensing images and the interferometry synthetic aperture radar, which show that the slope started to deform as early as 2011 (Fan et al., 2019).(Tian et al., 2020).The black line represents the mean value of annual rainfall.

Topographic Data
On 12 October, emergency teams performed uncrewed aerial vehicles (UAVs) aerial photogrammetry after the landslide.The KY-A001 UAV carried an HDR camera (with a frame of 23.4 × 15.6 mm and 5,210 pixels) to take the images.The Skyline Photomesh software (https://www.skylinesoft.com/photomesh/)was used to implement aerial triangulation and produce mosaicked images and dense point cloud data.We obtained a digital elevation model (DEM) with 10 m grid spacing before the landslide from the Sichuan Bureau of Surveying Mapping and Geoinformation (SBSMG).A 1.0 m DEM after the landslide was derived from the UAV-based photogrammetry (detailed processes described by Ouyang et al. (2019)).We down-sampled the 1.0 m DEM to 10 m grid spacing to obtain uniform DEM resolution before and after the landslide.With the two DEMs, we were able to determine the sliding surface and the 3D shape of the sliding body at its original location.Because the sliding material completely vacate the source zone, the sliding surface was constructed by taking the lower elevation from these two DEMs at each grid point.The shape of the sliding body was computed by comparing the sliding surface with the DEM before the landslide.
The field investigation showed that the landslide has a longitudinal length of about 1,500 m, an average width of about 600 m, and an elevation difference of about 850 m between the leading and trailing edges (Figure 1c).The sliding material from the source area has a total volume of about 24 million m 3 and traveled on an average slope angle of 24° along the direction of N95°E down to the valley, completely blocking the Jinsha River (Figure 1c).Finally, a landslide dam with about 2,000 m in length, 450 m in width, and a maximum thickness of 94 m formed and impounded a lake of about 290 million m 3 .The high-speed debris flowed horizontally for 1,200 m and then ran up about 160 m at the left bank of the Jinsha River (Figure 1c).

Seismic Data
Several seismic stations recorded the seismic signals from the 2018 Baige landslide.Zhang et al. (2019) selected five of these stations within distances between 108 and 255 km and filtered the significant long-period signals in the range of interest (30-140 s) using a jackknife method.The jackknife analysis approach is a fairly common practice to assess the uncertainty of an inversion (Moretti et al., 2015;Toney & Allstadt, 2021).Finally, the closest station was chosen to perform waveform inversion, and the three components of the impact forces on the surface were obtained.Here, we briefly introduce the basic theory of force inversion from seismic data.
A landslide event generally consists of an acceleration and then a deceleration phase.The shearing and loading applied by the landslide to the bed surface generate seismic waves.Based on the assumption of a single-force mechanism within the long-wavelength limit (Allstadt, 2013;Moretti et al., 2012;Yamada et al., 2013;Zhang et al., 2019;Zhao et al., 2015Zhao et al., , 2020)), the displacement (d) at a given seismic station for a seismic source can be calculated by convolving the force-time function (F) with the Green's function (G).In this case, d, F, and G are functions of time.The force-time function can be solved by deconvolving Green's function and long-period signals as follows: where G* represents the Green's function convolution matrix, G* T is the transpose of G*, α is the damping parameter chosen as giving the optimal trade-off between model norm and variance reduction of the seismic data.
For more details about this method, we refer to Allstadt (2013) and W. Li, Chen, et al. (2019).
Based on the assumption of a landslide as a rigid block sliding down a slope, the magnitude of the slope-parallel force results from the driving force of gravity opposed by the frictional force (Ekström & Stark, 2013;Moretti et al., 2015;Zhao et al., 2015).Thus, the apparent basal friction coefficient (μ) can be estimated from the mass and trajectory of sliding materials (Allstadt, 2013;Li, Chen, et al., 2019): where θ is the slope angle, g is the Earth's gravitational acceleration.f is the slope-parallel force exerted on the sliding block, which can be obtained from F. For more detail about seismic data inversion, please refer to Text S1 in Supporting Information S1.

Numerical Models
Numerical methods can simulate the landslide sliding process (including both the runout and deposit stages) and provide detailed kinematic information, including flow depth, velocities, and impact force.In contrast to depth-integrated methods, which cannot simulate the complex deformation pattern, the zones of shear concentration, and the momentum conversion mechanism inside the sliding body, the 3D stress-strain model can track the internal strain and stress and calculate the anisotropic stress distributions at each time step.This capability can significantly impact the extent and shape of the deposition.In this study, we will employ a 3D numerical method, smoothed particle hydrodynamics (SPH), with two constitutive models for the Baige landslide.We will briefly introduce the model and method in this section.

Smoothed Particle Hydrodynamics (SPH)
SPH is a meshless continuum-based method in which the mass and momentum conservation equations are solved based on a set of Lagrangian particles (Bui et al., 2008).Each particle carries physical variables, such as mass, density, and stress, and moves with a certain velocity.Therefore, the landslide is simulated in space and time by tracking each Lagrangian particle and its attributes in the computation domain.The simulations were performed using the 3D GPU-accelerated SPH solver LOQUAT (Peng et al., 2019), which assumes the sliding mass is homogeneous.In LOQUAT, the discretized governing equations for this model (please refer to the detailed deduction in Text S2 in Supporting Information S1) are: where d(⋅)/dt denotes the time derivative.ρ, m v, σ, and g are the density, mass, velocity, Cauchy stress tensor, and body force, respectively.i and j represent the so-called centering and neighboring particles.∇ is the gradient operator.n is the number of particles in the computation domain of the kernel function.W ij = W(x i − x j , h) is a smoothing kernel function, which depends on the distance between the two particles x i and x j , and the smoothing length h equals to 1.35 times particle size.The Wendland C2 kernel function (Wendland, 1995) has been employed in this study. Π = Π   + Π   is the stabilization term consisting of the artificial viscosity  Π   and artificial pressure  Π   , which are employed to damp out non-physical numerical oscillations and mitigate the tensile instability.I is an identity matrix.The second term on the right-hand side of Equation 3 is a numerical diffusion term used to improve the accuracy of pressure computation, where δ is a constant usually taken as 0.1 (Marrone et al., 2011), and ψ ij is the diffusion term.Equations 3 and 4 are ordinary differential equations that are solved using the second-order predictor-corrector explicit time integration method (Peng et al., 2019).In our SPH modeling, the non-slip rigid boundary was implemented using boundary particles (occupying a width of 2h) carrying the same initial spacing, mass, and density as fluid particles.These boundary particles participate in the SPH computation of the fluid particles like normal particles, however, the field variables are not computed by evolving the governing equations but extrapolated from the fluid domain.The boundary force f b is used to compare with the seismic inversion results and can be calculated by where i is the material particles.Π ai is the artificial viscosity term which is calculated by substituting subscript a for j.Note that Equations 3 and 4 do not include the computation of the stress tensor, which links the deformation and the mechanical response of the material (i.e., the constitutive relationships, which will be introduced in the following part).

Constitutive Models
We tested two different ways of computing the stress tensor using two constitutive models in the SPH simulations: the elastoplastic Drucker-Prager (DP) model from soil mechanics (Bui et al., 2008;Peng et al., 2019) as well as its counterpart in fluid mechanics, that is, the non-Newtonian rheological Drucker-Prager (RDP) model (Huang et al., 2012;Peng et al., 2022;Yang et al., 2020) (please refer to the detailed deduction in Text S2 in Supporting Information S1).
In 3D simulations, the Drucker-Prager (DP) model and the Mohr-Coulomb (MC) model are very similar, the only difference is that the yield surface of the MC model is a prism cone with shape edges, whereas the DP model is a perfect cone without any sharp edge.The DP model is a simplified version of the MC model, as it is more convenient to implement in numerical codes as it does not have sharp edges.As an elastoplastic model, the DP model includes the yield surface (f s ) and plastic potential (g s ): where J 2 is the second invariant defined as s: s/2 with s = σ − (I 1 /3)I, being the deviatoric stress tensor.
is the first invariant of the Cauchy stress tensor.ψ is the dilatancy angle controlling the volume change when soils are under shear deformation.k ϕ and k c are model parameters related to Mohr-Coulomb model parameter cohesion (c) and friction angle (ϕ) (Bui et al., 2008).The details on computing Cauchy stresses for the DP model can be found in Bui et al. (2008) and Zhan et al. (2019).
Different from the DP model, the RDP model assumes that the sliding material is a viscoplastic fluid with a yield strength (c + p tan ϕ).The material remains at rest while the stress state is below a threshold yield stress.Once the applied shear stress reaches the threshold yield stress, the mass starts flowing with the apparent viscosity (Becker & Idelsohn, 2016;Peng et al., 2021): where p is pressure,   is the scalar measurement of the strain rate tensor defined as  √ ̇γ ∶ ̇γ∕2 .The strain rate   γ = ∇ + (∇) T is two times of the symmetric part of the velocity gradient.m is a regularization parameter used to avoid the numerical singularity in the quasi-static regions when   has a very small value.A large regularization parameter leads to an extremely large apparent viscosity in the quasi-static flow regime that will severely reduce the numerical efficiency.Several previous test cases (S.Li et al., 2020;Peng et al., 2022) indicate that m = 100 s −1 can correctly simulate landslides.Therefore, in this study, we will set m = 100 s −1 with the RDP model for the Baige landslide.In addition, the pressure (p) can be obtained through an equation of state (EOS) based on the weakly compressible assumption (Bui et al., 2008;Huang & Dai, 2014;S. Li et al., 2020;Peng et al., 2019).
Both the DP and RDP models have two model parameters for describing physical properties related to stress: c and ϕ.Both models are based on the Coulomb friction law, assuming that the mass starts to flow when the shear stress exceeds the maximum shear strength provided by friction and cohesion.The most significant difference between the two models is that the DP model allows dilatation during shear deformation, whereas the RDP model assumes that the material is almost incompressible.It should be noted that both DP and RDP models describe the stress-strain relations of geomaterials in three dimensions.Therefore, both models consider the interaction between bedrock and sliding mass and the deformation and mechanical response of the sliding body, which can capture complex deformation and flow patterns in the sliding body.

Model Implementation
In this work, the formulations and models are implemented in the SPH solver LOQUAT developed by Peng et al. (2019).Large domains and volumes are common in landslide events, which often lead to a large number of SPH particles in 3D simulations.For instance, we used about 3 million particles for the Baige landslide with a 3.3 m diameter in the numerical domain.Simulations with large numbers of particles cannot be easily achieved using standard GPUs or parallel computing.The hardware acceleration using GPUs is implemented based on the Compute Unified Device Architecture (CUDA) framework, which is a programming and computing platform released by NVIDIA to enable scientific high-performance computing on GPUs.With GPU acceleration, our SPH solver LOQUAT can achieve more than a thousand times of acceleration compared to normal CPU-based serial computing (Chen et al., 2020;Peng et al., 2019).We can simulate high-resolution 3D cases for large-scale landslides on common desktops based on this.The simulation of the Baige landslide with around 3 million particles and 100 s of physical time can be finished in one to 2 hr, depending on the frequency of saving output files (i.e., see Text S3 in Supporting Information S1).The run times for the GPU-accelerated 3D SPH model are efficient, which is similar to those for a similar scale of the depth-averaged model using serial computing.However, this method can provide more detailed simulations in three dimensions, which makes it a promising alternative to serial-computing depth-integration techniques.
For the Baige landslide, the bedrock and the sliding materials are represented using SPH particles placed in a regular lattice.The particles were generated based on the topography data before and after the landslide, which is available in DEMs.The DEMs of the Baige landslides have a resolution of 10 m.We first interpolated DEMs to a 3.3-m resolution.Then the boundary particles were generated by converting the grid points in the DEMs into SPH particles and assigning the corresponding coordinates and properties, such as mass, volume, density, and stress.We added two additional layers of particles below the original layer to fulfill the requirement for the boundary of SPH simulations (Peng et al., 2019).The sliding body is then filled with particles into a regular lattice.Figure 3 shows how we generate the simulation domain for the Baige landslide.
Two major model parameters (cohesion and friction angle) for the rheological relations were tested by the trial-and-error method to calibrate the models.We optimized the parameters combination by a grid search similar to Yamada et al. (2018) to make sure the simulation results can be the best fit with the field observation data.This optimization method has been used successfully for the analysis of several landslides and debris flows (Coussot, 2017;Dai et al., 2014;S. Li et al., 2020;Yamada et al., 2018).We run a two-dimensional grid search for both the DP model and the RDP model with the following parameter space: c = 0.1, 1, 5, 10, 15, 20, 25, 30 kPa and tanϕ between 0.11 and 0.30 with a 0.01 or 0.02 increment (see Table S1 in Supporting Information S1).The choice of parameter ranges corresponds to the parameter values commonly used in large-scale landslides (Peng et al., 2022;Yamada et al., 2018), which can be expected to include the optimal parameter set.

Model Evaluation Indicator
To quantitatively evaluate model performance, we define some dimensional and non-dimensional parameters.For the topographic difference between numerical simulation and field observation, the following three parameters were used.The first parameter is the relative error of deposition area, η, defined as: where A s and A o are the deposition area of simulations and field observations, respectively.For an optimal model, η is equal to zero.
To overcome the weakness of the relative error of the deposition area in not separating false positives from false negatives, we used the point-to-point elevation metric (i.e., the mean squared error [MSE] and the normalized squared error [NSE]) and the time-force function metric to compensate for this deficiency.The MSE of the topographic thickness is defined as: where z s,k and z o,k are the thickness of the deposit at point k from the simulations and the field observations, respectively.n is the total number of grid points.The normalized mean square error, that is, the Nash-Sutcliffe efficiency (NSE) is defined as: where μ o is the mean value of the topographic thickness of the field observations.
The two metrics (MSE and NSE) evaluate the same information and are relative, for example, the minimum MSE matches the maximum NSE.However, MSE is dimensional, while NSE is dimensionless.So the authors left both metrics for readers to make it easier to read.
Unlike the field observation for topographic differences based on spatially distributed data, the seismic waveform inverted time-force function is a time-series data set.Here, the results of impact force from seismic wave inversion in the period range 30-140 s conducted by Zhang et al. (2019) are used to compare the forces calculated by numerical models.To evaluate the degree of agreement, we calculate the correlation coefficient (CC) between the two time-force curves (Moretti et al., 2015), as well as their normalized residual (NR) (Yamada et al., 2018) as follows: where f t and s t are the two forces at time t from the seismic wave inversion and the simulations, respectively.  and   are the mean values of the forces from the two sources.N is the total number of time points.The possible ranges of CC and NR are  [−1.0, 1.0] and  [0.0, +∞] , respectively.A perfect matching is implied by CC = 1.0 and NR = 0.0.But the range of η of the DP model (−51%-57%) is larger than the range of η of the RDP model (−54%-22%).We noticed a similar contour distribution across the parameter space of both DP and RDP models.This indicates that the minimum η can be achieved using different cohesion and friction coefficients combinations.In other words, the relative error between the simulations and observations could be zero, although with a differently shaped deposit.On the other hand, both the MSE and NSE show a nonlinear relationship with the cohesion and friction coefficient.Moreover, the parameter space within minimal MSE also has the maximal NSE.Based on the deposit fitness metrics, the ranges of parameters for DP and RDP models are about c = 13-25 kPa, tanϕ = 0.17-0.21and c = 1-15 kPa, tanϕ = 0.15-0.20,respectively.The point-wise comparison of topography differences helps find the optimal ranges of parameters for the two models.The gradient of the contour lines along the friction angle is larger than the gradients along cohesion for both the DP and RDP models (Figures 4a1 and 4a2).The gradients of MSE and NSE in the simulations with the DP model are more than 10 times larger than those with the RDP model when cohesion is less than 10 kPa (Figures 4b1, 4b2, 4c1, and 4c2).However, the large gradients are not observed in the RDP model (Figures 4b2 and 4c2).Generally, the MSE and NSE variation ranges from the RDP model are smaller than the ranges of MSE and NSE of the DP model.

Comparisons Between Modeled and Observed Topography Difference
Within the parameter space from the optimal ranges, we compared the best simulation (c = 15 kPa and tanϕ = 0.19 for DP model and c = 15 kPa and tanϕ = 0.15 for RDP model) and the observed sediment spatial distribution, the topography difference, and the distribution of elevation changes at each grid cell before and after the event (Figure 5).The elevation changes in the deposition area were calculated from the pre-failure DEM (2017) and post-failure DEM (12 October 2018) (Figure 5a).The deposition area was oriented parallel to the river channel, and the maximum deposit thickness was 94 m.The debris even propagated to the right bank of the Jinsha River and scoured the toe of the slope (see the dark gray color in Figure 5a).Generally, the horizontal extent (area and runout) and the deposit thickness obtained by the numerical simulations agree with the field investigation (Figures 5a-c).However, some notable differences appear in the landslide deposit thickness.For example, some SPH particles stayed on the slope after sliding, resulting in region-thin deposits of 1-2 layers (5-6 m) of SPH particles (Figures 5b and 5c).The simulated maximum thickness and its adjacent range do not agree well The deposition area includes the dammed lake and the sliding surge impacted area on the opposite slope but excludes the downstream area.Note, the artificial bottom boundary (replaced by river surface elevation due to the lack of underwater terrain) at a regional dammed lake may result in some differences in calculating regional deposit thickness.
with the field observations.The frequencies of the numerical results are all shifted to the right side of the field observations (Figure 5d), in other words, the simulations tend to over-predict the local deposit thickness, for instance, the sliding slope and the area downstream of the blocked river.The overestimated region is on the left bank of the river (see the wave-affected area in Figure 1b and the right-bottom deep blue range in Figure 5a), while the underestimated region is at the dammed reservoir (see the red dashed lines in Figure 5a).The detailed point-wise comparison of deposition thickness also shows overestimation at the slope and the downstream part of deposition but underestimation at the upstream part of the deposition (Figures 5e and 5f).

Comparisons Between Modeled and Inverted Impact Force History
Figures 6 and 7 show the CC and NR of the DP model and the RDP model in the two-dimensional parameter space.The gradient of the contour lines for CC and NR are almost parallel to the x-axis, which implies a low sensitivity to cohesion (Figures 6 and 7).The vertical component CC is closer to 1.0 (NR is closer to 0.0 accordingly) compared with that for the east and north components.Moreover, both CC and NR of the east and north components decrease monotonically with increasing cohesion and friction angle but have non-linear relationships with respect to cohesion and friction coefficient in the vertical component.The highest value of CC or the smallest values of NR in the vertical component are mainly found around c = 5-20 kPa and tanϕ = 0.16-0.20 for both the DP and the RDP models.
Figure 8 shows the comparison between forces in the three different components from the simulations with the optimal parameters of the DP model and RDP model and the inverted time-force curve from the seismic signals.The inverted force-time function within 80 s duration of the landslide movement is shown in Figure 8 (see the green lines).The force at the vertical component showed a downdip direction within t = 0-20 s, which represents the acceleration stage of the landslide source that suddenly slides.Subsequently, a deceleration stage (t = 20-40 s) happened when the slide mass reaches the river channel; thus, the force reversed direction compared to the preceding stage.Then, the energy of the landslide decreased to a relatively small level, accompanied by another acceleration process (t = 40-57 s) and deceleration process (t = 57-80 s), which are generated by the landslide impacting the opposite hillside when the flowing mass climbs up the left slope and then spreads far away in the up-and downstream directions associated with the depositing stage of the landslide.The simulated amplitude and phase of the forces correspond closely to the inverted force-time history.For example, the simulated first two waves (t = 12-15 s and t = 29-31 s) compare well to the inverted results on the vertical and east components.The amplitude of the first peak of the east component from the simulation is higher than the inverted force.Also, some mismatches exist from approx.40-80 s, which may be caused by the landslide impacting and eroding on the opposite hillslope and the interaction between landslide and river flow that we do not account for in the simulations.Meanwhile, there are phase shifts (approx.15 s) between the numerical and seismic inverted results on the east and north components.In other words, the initiation of the simulated forces at the east and north components lag behind that of the seismic inversion.The discrepancies may stem from the single-force mechanism that does not include additional aseismic transport of material (Hibert et al., 2014).In addition, we do not take into account that the rockslide mass disintegrates and fragments in the initial stages of motion.

Model Performance
Our models generally represent the landslide event regarding the deposition area, the deposition height, and the runout distance (Movie S1 Supporting Information S1).Differences between model results and observations are mainly in the spatial distribution of sediment thickness.In other words, the comparison of deposit morphology measured in the final static deposition does not allow us to determine the best model parameters.The comparison of point-to-point deposition, as well as the dynamic time-force function, help us to compensate for this deficiency.That is to say, a comprehensive and effective approach to numerical model calibration is to compare simulation results to both static field observations and dynamic seismic signals.It should be mentioned that when using seismic inversion data to constrain runout models, it is best to use multiple stations and to assess the inversion uncertainty when interpreting results.
The time-force sequence from our simulations is consistent with the inverted force from seismic signals in amplitude and phase.The simulated and inverted force in the vertical component match better compared to the other two components (Figures 6-8).This is because we adjusted the representative force of the slide mass to ensure the vertical impact scales fitted to the seismic inversion, and excluded the agreement in the horizontal direction.The discrepancies in the three components indicate complicated spatial behavior, for instance, in slide mass fragmentation, acceleration, deceleration, impacting the opposite hillslope and rebound, blocking the river and dam formation, and so on.Compared with the deficiency of ignoring the variation of vertical shear stress in the two-dimensional model, the 3D stress-strain relationship SPH formulation can capture velocity, strain, and stress variation inside the sliding body at each time step because it considers the friction interaction between the sliding body and the bedrock and within the interior of the sliding material itself.For instance, the internal anisotropic stress distributions are tracked, that is, if the flow converges in one principal direction and diverges in the other.This capability can have a significant impact on the extent and shape of the deposition.Therefore, the complex deformation pattern (e.g., forward motion, lateral spread, and vertical runup) was well simulated.Other benefits of 3D simulations are, for example, (a) with 3D simulations, the dynamic fragmentation, the zones of shear concentration, and the momentum conversion mechanism inside the sliding body can be more accurately simulated.For example, a landslide may occur progressively, with some part of the sliding body failing first, and then the failure propagates backward and downwards till the bedrock is reached.With 3D simulation, this complex failure behavior can be reproduced and analyzed; (b) 3D continuum-based simulations serve as the base of more complex analysis, such as considering the influence of huge rocks and tree trunks in a continuous-discrete numerical framework; (c) last but not least, with 3D simulations the interaction between landslides and structures can be more accurately analyzed.However, we are fully aware that our current work does not go so deep to demonstrate the benefits of employing 3D simulations comprehensively.In addition, the high-efficient GPU acceleration is specifically beneficial in the varied and complex topography of real-world terrain (Peng et al., 2022).
Although both the DP model and the RDP model can reproduce the landslide processes with reasonable parameters, some differences are observed in the deposition patterns between the two constitutive models.To investigate these differences, we plotted the velocity distributions (Figure 9) at the cross-section A-A (see the dashed line in Figure 1c) during the slide progression.The velocity magnitude of the DP model is almost larger than that of the RDP model, for example, the maximum surface velocity of the DP model reaches about 75 m/s and 50 m/s while RDP only has approx.65 and 30 m/s at t = 30 and 50 s, respectively.This means that the impact between the already deposited material and the flowing material is stronger in the simulation with the DP model, which is the reason that a larger final deposition area is obtained with the DP model.Specifically, a significant climb-up occurs in the DP model, while the profiles of the RDP sliding are relatively flat at t = 50 s.The reason may be that the numerical volumetric expansion in the DP model (Yang et al., 2020) give rise to faster material movement and larger expansion than RDP model.This non-physical volumetric expansion does not change the friction angle but is more like a numerical deficiency in large deformation analysis using SPH or other particle-based numerical methods.Moreover, the time-series velocity and stress at the cross-sectional of sections A-A and B-B (see Figure 1c) are plotted in Figure 10.The velocity distributions from the simulations using both models follow approximately the Bagnold velocity profile (Hill et al., 2022;MiDi, 2004).The maximal velocities first increase and then decrease, with a peak velocity over 60 m/s.The cross-section velocity from the DP model is also larger than that from the RDP model for most of the flow process, which is consistent with our analysis that the DP model has a larger non-physical volumetric increase.As for the stress, results from the two simulations show mostly linear increases with respect to the material depth, except for the initial phase.The stress fluctuation of the DP model is smaller than that of RDP, which is the result of their ways of computing the stress tensor and the treatment of plasticity is different.The DP model is a pure elastoplastic model, which derives the stress tensor over time using a stress-strain relationship that relates the stress increment to the strain increment.Instead, the RDP model is a non-Newtonian model, or a so-called viscoplastic model, which evolves the stress tensor from instantaneous density and strain rate.The DP model can consider the elastic regime, while the RDP model assumes that the materials are always under plastic flow.Therefore, there are strong stress fluctuations in the profile produced by the RDP model.In addition, the DP model allows the sliding volume to dilate, but the RDP model does not.In the DP model, the volume can increase depending on the stress status during the flow process while the mechanical properties described by the model remain the same.In the RDP model, the material is Figure 6.Results of numerical simulations for two calibration parameters (Friction angle, tan ϕ, and cohesion, c) as contour maps based on seismic signal inversion results using correlation coefficient: (a1) and (a2) are the correlation coefficient (CC vertical ) between vertical impact force from seismic wave inversion and our simulations.(b1) and (b2) are the correlation coefficient (CC east ) between impact force in the east direction from seismic wave inversion and our simulations.(c1) and (c2) are the correlation coefficient (CC north ) between impact force in the north direction from seismic wave inversion and our simulations.The results for the DP model are in the left column, and the results for the RDP model are in the right column.The black circular on each panel of the figure marked the best-fit parameters.assumed to be weakly compressible (usually by less than 1%), which ensures that the volume of the sliding material is generally unchanged in the whole sliding process.In general, the overall lateral spreading of the final deposit predicted by the RDP model is much smaller compared to the DP model (see Figure 1 in Supporting Information S1), which is consequential of the different stress treatments and fundamental assumptions used in the two models.

Model Parameter
The motion and deposition of the sliding mass are mainly controlled by complex mechanics reflected in the rheological features.We calibrated the model parameters (cohesion and friction angle) by comparing modeled results with field observations and seismic signals.In the range of geomaterials cohesion (0.1-30 kPa), the simulated force-time function (e.g., CC vertical changes from 0.01 to 0.14 for the RDP model) and deposition area (e.g., η changes from 8% to 18% for RDP model) are relatively limited in the simulation.The simulation results indicate that both the deposit extension (Figure 4) and time-force sequence (Figures 6 and 7) are more sensitive to the friction angle than the cohesion.In other words, in the large-scale Baige landslides, the major contribution of the yield stress is friction rather than cohesion.This is because, in large-scale landslides, the yield stress is mainly dependent on friction caused by the high normal stress due to the large thickness of sliding mass.Therefore, we can also infer that cohesion generally controls the final extension and cessation of the edge of deposits because of the small material height close to the margin where the friction influence is too small to eclipse the effect of cohesion (Peng et al., 2022).
The aforementioned results are like the difference between small-scale and large-scale problems.In small-scale problems, cohesion is very important and largely determines the final deposition (Bui et al., 2008;Peng et al., 2022).In contrast, in large-scale problems, cohesion plays a lesser role.The potential reason for this is that cohesion depends on the landslide scale, whereas the frictional force is proportional to the normal contact force.Based on this, we can speculate that friction is more critical than cohesion in large-scale deep-seated landslides.

Shear Strength
For the Baige landslide, we assume that the rock and soil behave as Mohr-Coulomb materials.Therefore, the main factors contributing to sliding are Mohr-Coulomb shear strength parameters (cohesion c and friction angle tanϕ).Comparing the simulated and field-observed results allows us to infer the shear strength properties during the Baige Landslide.The optimal range of friction values is relatively narrow (0.19-0.21 for the DP model, 0.12-0.15for the RDP model).These friction values are much smaller than the depth-averaged method (0.30-0.84,McDougall & Hungr, 2004;Ouyang et al., 2019;Pastor et al., 2009Pastor et al., , 2014)), the repose angle from field investigations (0.36-0.78,Tofani et al. (2017)), and laboratory test results (0.27-0.50, De Blasio et al., 2011;Kaitna et al., 2007;Parsons et al., 2001), which are commonly used in small-scale landslide modeling (Bui et al., 2008;Peng et al., 2019).However, the modeled friction coefficient is in agreement with those inferred from large-scale landslides, in which values are as low as 0.10-0.27(Aaron, 2017;Duncan et al., 2014;Keefer & Johnson, 1983;Lucas et al., 2014).One of the potential reasons is the velocity-dependent friction weakening that is hypothesized to happen during the sliding process of the giant landslides (Lucas et al., 2014).This indicates

Model Limitation and Future Works
Landslides are complex multiphase flows, in which possibly multiple mechanisms are active at different stages.Thus, it may be thought that an ideal model takes all phases and their interactions into account.Indeed, recently, models of geophysical flows have developed to incorporate an increasing number of phases and components (e.g., Hungr & McDougall, 2009;Iverson & George, 2014;Leonardi et al., 2014;Mergili et al., 2017;Pudasaini, 2012;von Boetticher et al., 2017).However, these representations are far from the complex reality of nature due to all the simplifications, especially for ignoring the relative motion between phases.Given the complex interactions between different parts of a sliding mass and that most of the model parameters cannot be obtained in the field, no universal numerical model can perfectly simulate all landslides.In this study, two simple homogeneous models were adopted as a first approximation to simulate the Baige landslide.Our models can give a relatively clear picture of the 3D dynamic processes for the landslide that resemble actual events.Yet, some limitations to our model still remain.First, the rheological relationship is based on a single-phase homogeneous assumption.However, large boulders and rock fragmentation effects have not been modeled and may influence the motion.Second, the model does not consider soil erosion on slopes and ignores the impact of river flow after the landslide, which can significantly change the sediment distribution.Third, the friction reduction (Lucas et al., 2014) at the sliding surface during the sliding process is not considered, which may underestimate the landslide mobility.In actual landslide processes, the equivalent friction usually decreases due to the pore fluid pressure under the high-speed shearing condition (Iverson et al., 1997;Schulz et al., 2009).Especially when the landslide mass reaches the bottom, the high pore fluid pressure due to the mass compression and interaction with the water will significantly change the deposition pattern.Our model does not consider the stress changes inside soils caused by the pore-water pressure changes and the progressive cracks during the initiation period.Unfortunately, we do not have a very clear understanding of how pore pressure changes during landslide processes.Therefore, advanced-scale experiments are needed to explore this relationship.Last, we only considered the pre-landslide surface elevation on the left of the river bank and assumed the post-landslide sediment surface was the same as the dammed lake upstream.These artificially increased the basal elevation at the left river bank and overestimated deposit thickness at the reservoir.

Conclusion
In this study, we mainly focus on implementing numerical simulations of a landslide.By comparing the model results of two constitutive models with both static final deposition data and dynamic seismic signals, we checked the sensitivities of the models and involved parameters.The main findings are as follows.
1.A continuum-based mesh-less method, the 3D GPU-acceleration SPH method has clear advantages in high resolution and low-cost computation to simulate large-scale natural landslides in complex terrain.Both the soil mechanic model (DP) and the non-Newtonian fluid model (RDP) can reproduce key dynamic processes and deposition morphology.Although some differences exist between the two models, their overall performance is comparable.2. The calibration indicators such as deposit area, runout, and maximum thickness are insufficient to judge the accuracy of models.A point-wise comparison of each elevation can minimize the error to give a more confident model.The conventional calibration-based method (numerical-and field-obtained final static data) can only estimate the bulk characteristic of landslides but provide limited information on the dynamic process.The seismic signals can provide important timing information complementary to the post-event deposit information.By utilizing both the static field observations and dynamic seismic waves, and multi-criteria calibration indicators, numerical models, can be better calibrated.3. The Mohr-Coulomb shear strength parameters (friction and cohesion) are crucial factors in controlling landslide motion and deposition.The friction coefficient has a larger effect on the sliding dynamics and final deposition compared with cohesion for deep-seated landslides.The friction coefficient obtained from 3D SPH modeling is much less than the depth-average method and laboratory tests, indicating that internal stress-strain contributes greatly to the material movement and there are unrevealed mechanisms leading to friction reduction.

Figure 1 .
Figure 1.The study site is located at the Baige village, the border between Jiangda county and Baiyu county, China: (a) Location of the study area in the southeastern part of the Qinghai-Tibet Plateau; (b) Zoom in the location of the Baige landslide, (c) Top view photographs of the landslide (image taken by a UAV on 12 October 2018 from Ouyang et al. (2019)) to show the source area, accelerated area, deposition area (landslide dam) and wave-affected area for this landslide.

Figure 2 .
Figure 2. Precipitation data of the study area from the Dege Country Meteorological Station: (a) Daily and accumulated rainfall from 2016 before the landslide.The red cross represents the date when the Baige landslide happened.(b) The annual rainfall from 1957 to 2018(Tian et al., 2020).The black line represents the mean value of annual rainfall.

Figure 3 .
Figure 3. Actual topographic difference during landslides and the smoothed particle hydrodynamics (SPH) particle generation procedure: (a) Remote sensing figure shows the landslide scare after the Baige landslide.The landslide body is estimated based on DEM difference before and after the landslide (b).We generate SPH particles based on the digital elevation model differences (c) and represent the sliding body (d).

Figure 4
Figure4shows η, MSE, and NSE of the DP model and the RDP model in the two-dimensional parameter space.Generally, with the increasing cohesion and friction angle, η gradually changes from positive to negative values.

Figure 4 .
Figure 4. Results of numerical simulations for two fitting parameters (Friction angle, tan ϕ, and cohesion, C) as contour maps based on deposition pattern (a1) and (a2) are the relative errors of deposition area (η) with a range of −50%-50%.(b1) and (b2) are the mean squared errors (MSE) with a range of 0-30 m (c1) and (c2) are Nash-Sutcliffe efficiency (NSE) for each model with a range of 0-1.The results for the DP model are in the left column, and the results for the RDP model are in the right column.The black circular on each panel of figure marked the best-fit parameters.

Figure 5 .
Figure 5.The deposit distribution of the Baige landslide.(a) Topography difference from field investigation and simulation results from (b) Drucker-Prager (DP) model (c = 15 kPa and tanϕ = 0.19) and (c) rheological Drucker-Prager (RDP) model (c = 15 kPa and tanϕ = 0.15).The black line shows the extent of the landslide.The white arrows indicate the flow direction in the river.The red dashed line is the boundary of the dammed lake after the landslide.The red line shows the landslide mass flow along the river during the sliding.The frequency and difference of the topographic changes between the numerical simulations and the field observation are plotted in (d), from the DP model (e) and RDP model (f).The deposition area includes the dammed lake and the sliding surge impacted area on the opposite slope but excludes the downstream area.Note, the artificial bottom boundary (replaced by river surface elevation due to the lack of underwater terrain) at a regional dammed lake may result in some differences in calculating regional deposit thickness.

Figure 7 .
Figure7.Results of numerical simulations for two calibration parameters (Friction angle, tan ϕ, and cohesion, c) as contour maps based on seismic signal inversion results using normalized residual: (a1) and (a2) are the normalized residual (NR vertical ) between vertical impact force from seismic wave inversion and our simulations.(b1) and (b2) are the normalized residual (NR east ) between impact force in the east direction from seismic wave inversion and our simulations.(c1) and (c2) are the normalized residual (NR north ) between impact force in the north direction from seismic wave inversion and our simulations.The results for the Drucker-Prager (DP) model are in the left column, and the results for the rheological DP (RDP) model are in the right column.The black circular on each panel of the figure marked the bestfit parameters.

Figure 8 .
Figure 8. Impact force exerted by the landslide (a, b, and c represent the vertical, east, and north directions): the different lines represent the inverted force from seismic signals (green), the best-simulated force corresponding to the Drucker-Prager (DP) model (red) and rheological DP (RDP) model (blue).The start time of simulations is zero, within the period of seismic signals.The fitting parameters for the DP model are c = 15 kPa and tanϕ = 0.19, and for the RDP model are c = 15 kPa and tanϕ = 0.15.Calibration parameters for vertical component: DP, NR vertical = 0.276 and CC vertical = 0.852; RDP, NR vertical = 0.363 and CC vertical = 0.804.Calibration parameters for east component: DP, NR east = 1.056 and CC east = 0.706; RDP, NR east = 0.769 and CC east = 0.798.Calibration parameters for north component: DP, NR north = 1.064 and CC north = 0.332; RDP, NR north = 0.946 and CC north = 0.394.

Figure 9 .
Figure 9. Time-series velocity magnitude at section A-A (see Figure 1b, about 1,500 m) for Drucker-Prager (DP) model (left) and rheological DP (RDP) model (right) with the best fitting parameter sets.The sliding body in the DP model propagates faster than those from the RDP model as best fitting cases.Meanwhile, the RDP model has thicker sediment than the DP model at the slope but has a smaller runup at the opposite slope.At t = 50 s, the sliding body in the DP model shows a shock wave behavior.In contrast, the sliding body in the RDP model is smooth.

Figure 10 .
Figure 10.Profiles of time-series velocity (left) and stress (right) are extracted at the cross-sectional of sections A-A and B-B.