Mechanism of Co‐Seismic Deformation of the Slow‐Moving La Sorbella Landslide in Italy Revealed by MPM Analysis

Predicting the seismic behavior of slow‐moving landslides presents a significant challenge. Their yearly displacements of a few millimeters to several meters indicate that even before any earthquake they are balancing at the verge of instability. It is therefore not surprising that conventional analysis predicts large co‐seismic displacements even for light‐to‐moderate ground motions. In reality, however, negligibly small displacements have often been observed for such earthquakes, while strong ground motions can still lead to catastrophic failure. This discrepancy challenges both our understanding of the underlying mechanisms and the reliability of conventional analysis, resulting in significant uncertainties for risk assessment. Progress in this field has been further hindered by incompleteness of the available monitoring data and a lack of reliable numerical models capable of dealing simultaneously with small and large deformations under seismic loading. In this article, recently published monitoring data on La Sorbella landslide in Italy and the latest developments in seismic material point method analysis are used to gain insights into the mechanisms controlling the co‐seismic behavior of slow‐moving landslides. This forms the basis for a subsequent investigation of the potential landslide behavior during strong motions, using various material models. The study demonstrates that the co‐seismic behavior of La Sorbella slide is strongly influenced by the interplay between geometrical and viscous effects, as well as potential softening in the shear zone and the surrounding soil mass. The proposed approach paves a way toward quantification of these effects in risk assessment for slow‐moving landslides.


Modeling the Seismic Response of Slopes and Landslides
Numerous techniques have been developed in the past for assessing the stability of slopes (Jibson, 2011), ranging from the pseudo-static limiting equilibrium approach to Newmark's sliding block analysis (Newmark, 1965), and to stress-deformation analyses (Hashash & Groholski, 2010;Kramer, 1996). Whereas the former allows the stability of the slope to be estimated during the seismic event, Newmark's sliding block analysis provides a method to assess the co-seismic displacements of a rigid block. However, for real case geometries and non-linear soil behavior, this simplified approach cannot be relied on to evaluate earthquake-induced slope deformations. For rigorous stress-deformation analysis, the finite element method (FEM) is the most frequently used technique and has been successfully applied for seismic analyses (C. Wang et al., 2019;Stoecklin et al., 2020;Zheng et al., 2005). However, for mesh-based methods (e.g., FEM), special consideration has to be used when dealing with large deformations, owing to the problem of mesh distortion. To overcome this difficulty, several alternative numerical methods have been proposed, in particular the coupled Eulerian Lagrangian finite element method (CEL) and the material point method (MPM). A comprehensive review of different numerical approaches was provided by Soga et al. (2016). With the CEL a traditional Lagrangian phase, in which elements deform with the material, is followed by an Eulerian phase, during which elements with significant deformation are remeshed and material flow between neighboring elements is computed (Dassault Système Simulia, 2019). The material is tracked and evolves through the Eulerian mesh based on its volume fraction, which entails a smearing of the state variables and therefore is less suitable for the distinct contrast given by shear zones in slow-moving landslides. The main disadvantage is that existing CEL implementations do not include appropriate seismic boundary conditions and hence cannot be applied for seismic response analyses. Traditional Lagrangian FEM does include appropriate boundary conditions (Nielsen, 2014) and could theoretically be used to model moderate displacements where large strains are concentrated in the shear zone. This requires that the shear zone is introduced using a contact formulation (Wriggers, 1995). However, this methodology comes with convergence and accuracy issues, which are exacerbated by the nonlinear geometry of a real shear zone. MPM, on the other hand, represents a suitable framework for the seismic response of slopes, because it combines the advantages of standard Lagrangian and large-strain Eulerian approaches, but until recently it has lacked a rigorous implementation of the seismic boundary conditions. This implementation was recently performed by Kohler et al. (2021), making the MPM framework applicable for investigation of co-seismic displacements of slowly moving landslides.

Material Point Method (MPM)
The MPM is a hybrid approach in which the material is represented as Lagrangian particles, while the equations of motion are solved on an Eulerian grid. The Eulerian solution procedure allows for the material to undergo large deformations, while the Lagrangian material representation provides a convenient way of tracking material properties and constitutive state variables (e.g., stresses and strains). MPM is a generalization of the Particle In Cell (PIC) and the Fluid Implicit Particle (FLIP) methods for solid mechanics and was first proposed by Sulsky et al. (1994). Nowadays, numerous variations of the original algorithm have been introduced, tailored for a range of different applications (Bardenhagen & Kober, 2004;Jiang et al., 2015;Steffen et al., 2008). In recent years, the MPM has been successfully applied to model seismic slope failure (Bhandari et al., 2016;Ering & Sivakumar Babu, 2020;Kohler et al., 2021) and landslides triggered by earthquakes (Alsardi & Yerro, 2021;He et al., 2019). Alvarado et al. (2019) used a hydro-and thermo-mechanically coupled MPM to model the non-seismic behavior of active landslides.

Goal and Structure of This Paper
As will be shown in the following section, a surprising feature of the co-seismic behavior of La Sorbella landslide is that its displacements recorded in 2016 during three significant earthquakes were extremely small (<1 mm). This is for an intrinsically unstable slope, which can move around 0.5 mm/day in the absence of earthquakes. There exist other cases where relatively small co-seismic displacements of active landslides have been reported (Al-Homoud & Tahtamoni, 2000;Bontemps et al., 2020;Lacroix et al., 2014). As will be shown in this paper, the conventional Newmark's sliding block analysis would predict, for a slow-moving landslide whose safety factor approaches unity, a displacement magnitude tending to infinity. To resolve this contradiction and to give some insight into the controlling factors, an effort has been made to develop a more rigorous modeling approach dealing with active landslides under seismic loading.
In the first part of this paper, La Sorbella landslide is introduced, and the necessary geological and geotechnical properties are discussed. Special focus is placed on its behavior during three recent earthquakes. This is followed by a section in which a modeling technique for assessing the co-seismic behavior of active landslides based on the MPM is presented. These landslides often feature distinct zones of intensive shearing which will be addressed by introducing the large-deformation constitutive model for shear zones. This approach offers the advantage of the same model being used both for lighter earthquakes resulting in small deformations and for stronger motions causing a non-linear soil response that leads to a catastrophic failure. For a mechanical and numerical validation of this technique, a benchmark model for infinite slope conditions is presented and compared to the conventional Newmark's sliding block analysis and a finite element approach. Applying this methodology to La Sorbella landslide for recorded moderate earthquakes will allow an understanding of the relative contribution of geometric effects, groundwater level and rate dependency of soil strength to the reduction in co-seismic displacements. Finally, the potential of the presented methodology will be shown by subjecting the model to stronger motions and investigating the influence of the shear zone behavior.

Geological Description
La Sorbella is an active deep-seated landslide in central Italy located near Valfabbrica (Figure 1a), which was already shown as active in the map by Guzzetti and Cardinali (1989). Monitoring of the landslide began in 2001, whereas intensive data acquisition through inclinometers started in 2014. All the necessary information and data concerning La Sorbella landslide have been adopted from Ferretti et al. (2019) and Ruggeri et al. (2020), where the monitoring campaign and measurement results are presented in detail. The landslide spans a slope of an average inclination of 8°. It has a length of 550 m and a maximal width at the toe of 600 m (Figure 1b). According to Ferretti et al. (2019), soil samples from the landslide mass consist of 70%-90% silt and clay and are classified as CL (clay of low plasticity) and ML (silt of low plasticity). Unfortunately, the available data from laboratory tests are rather limited. According to Ferretti et al. (2019), cohesion in the sliding layer is negligible, the peak friction angle ranges from 28° to 32° and the residual friction angle lies between 14° and 18° (Table 1). No tests were performed on the material from the shear zone, but the residual friction angle is assumed to be smaller than 14°, otherwise the slope would remain stable. This assumption can be justified by the fact that the material in the shear zone is likely to be weaker than in the sliding layer and has undergone very large deformations. A possible range of the residual strength in the shear zone will be back-calculated in Section 4.3 using a strength reduction analysis. The underlying bedrock in this region can be assigned to the Marnoso-Arenacea formation, which is characterized by alternating arenaceous and pelitic layers (Guerrera et al., 2012). Geophysical properties are not directly available for La Sorbella landslide and are assumed based on cross-and downhole seismic tests from a similar landslide in the Marnoso-Arenacea formation (Assefa et al., 2017). These assumptions are confirmed by the geophysical profile available for the seismic station IV-ATCC (Russo et al., 2022) which is located 6 km from the landslide and reaches the Marnoso-Arenacea formation.
Based on inclinometer measurements, whose locations are shown in Figure 1a, the currently active landslide can be clearly identified from a distinct shear zone of less than 1 m thickness, where almost all the deformation takes place (Ruggeri et al., 2020). The shear zone runs mostly at a depth between 20 and 36 m along the transition from the bedrock to the overlying landslide mass. A scarp at the top of the landslide and no significant deformations at the inclinometer AI8 at the toe indicate that the shear zone runs through the landslide mass to the surface at the lower and upper ends (Ruggeri et al., 2020). Moreover, the inclinometer measurements ( Figure 1c) suggest that the landslide moves largely as a rigid body with little internal deformation and at a very slow rate of 1.0-1.5 cm/year (Ruggeri et al., 2020). According to Ferretti et al. (2019), the groundwater level oscillates seasonally between 2 and 6 m below the surface, but the actual piezometer data are not available. In the study by Ruggeri et al. (2020), it is clearly shown that the landslide displacements are closely linked to rainfall events.

Seismic Behavior
In 2016, a long series of earthquakes occurred along an Apennine fault system located around 50 km southeast from La Sorbella landslide (Figure 1a). In this sequence, three events were characterized with a magnitude ≥ 5.5 (Table 2)

of 25
located less than 1 km from the landslide. Only limited information is available about the site characteristics of the station. The seismograph is installed in the basement of a building, and the site is classified based on the surface geology as ground type B (deposit of very dense soil at least several tens of meters in thickness and a gradual increase of mechanical properties with depth, v s.30 = 360-800 m/s) according to Eurocode 8 (Comité Européen de Normalisation, 2004). In the station report (Russo et al., 2022), the average shear wave velocity in the top 30 m is estimated as v s.30 = 550 m/s. The H/V spectral analysis shows only an indistinct peak at a frequency of 4.8 Hz for ambient vibrations according to SESAME (European Commission, 2004), and the seismic H/V analysis for 49 earthquake motions did not reveal any peak or eigenfrequency (Russo et al., 2022).
The displacements recorded by an automatic in-place inclinometer system located in borehole AI11 and the precipitation data from the Casanuova dam weather station (Servizio Idrografico Regione Umbria, 2016), located 2 km from the landslide, are presented for the days around the earthquake series in Figure 1c, which has been adapted from Ferretti et al. (2019). The co-seismic displacement corresponding to the recorded earthquakes can be clearly distinguished, as well as the fact that the landslide has been in an active state of slow movement during this period. Detailed analysis by Ruggeri et al. (2020) of the data from several years of measurements demonstrated that, except for these three earthquake events, a maximal displacement rate of 0.4 mm/day is reached only during periods of heavy precipitation over several days. Another interesting fact pointed out by Ruggeri et al. (2020) is that no change in sliding regime after the seismic events could be identified. In contrast, Lacroix et al. (2014) and Bontemps et al. (2020) reported that in the case of the Maca landslide in Peru, the co-seismic displacements are followed by three times greater post-seismic slip over the following 5-week period. The lack of post-seismic acceleration in the La Sorbella case can be explained by the presence of a distinct and well developed shear zone, where the clayey soil reached the residual state with shearing taking place at constant volume, without changes in pore water pressure (Skempton, 1985).

MPM Framework
The basis for modeling the seismic response of landslides is the MPM framework proposed by Kohler et al. (2021), where the implementation of appropriate boundary conditions for the seismic model is presented. The MPM code  (2015) and Stomakhin et al. (2013) and can be described in a four-step algorithm, shown in Figure 2.
The seismic response of slopes is modeled in two main steps (Kohler et al., 2021). First, the static stress field within the slope is computed using kinematic boundary conditions. In a second step, the actual seismic simulation is performed by applying the concept of a compliant base (Lysmer & Kuhlemeyer, 1969) and free-field columns (Wolf, 1989;Zienkiewicz et al., 1989). In order to model the existing landslide, the shear zone is pre-defined by assigning the corresponding material points to a different constitutive model (see Figure 2b). The thickness of the shear zone has to be chosen based on the span of the interpolation functions to allow the material points to move accordingly. For cubic B-splines, a thickness of at least double the grid size is recommended. Except for very small landslides, it is not reasonable to match the numerical shear band thickness to the in-situ zone of intensive shearing, as a very small grid size would be needed. The computational cost would be rather high, owing to the small stable time increment and the large number of material points in the case of a regular grid. However, the numerical thickness has to be taken into account for the constitutive model, and an appropriate scaling approach has to be applied. Therefore, the smeared crack approach introduced by Rots et al. (1985) is often applied in MPM (Kohler et al., 2021;Soga et al., 2016;Yerro, 2015).

Constitutive Model for the Shear Zone
For modeling a creeping landslide, an appropriate constitutive model for the shear zone needs to be implemented. Considering that MPM is a large-strain solution method, the choice of stress-and strain measures should be carefully assessed in terms of objectivity. The seismic behavior of active landslides must be simulated regardless of the magnitude of displacements since they might vary significantly depending on the initial conditions, material parameters and input motion. It has been shown that the straightforward extension of small-strain models to finite strains based on rate form equations leads to various inconsistencies, such as spurious stress oscillations or improper energy dissipation (Bažant et al., 2012;Perić et al., 1992;Simo & Pister, 1984). Therefore, it is advantageous to use the multiplicative decomposition of the deformation gradient into an elastic and plastic contribution, as introduced by Lee (1969) = where and are the elastic and plastic deformation gradients. A convenient expression for describing the elastic deformation is provided by the Eulerian logarithmic strain or Hencky strain (Hencky, 1928) where denotes the elastic left Cauchy-Green strain tensor. The elastic constitutive relation is given by the energy density as with the Lamé constants and. The Kirchhoff stress can be defined analogously to the small-strain approach as where represents the identity tensor. For the framework presented here, a viscoplastic material based on the consistency model (W. M. Wang et al., 1997) is applied, where, in contrast to over-stress models, a rate-dependent yield surface is introduced. Based on the constitutive model proposed by Wedage et al. (1998), a Drucker-Prager yield surface matched to a Mohr-Coulomb surface in plane strain for flow at constant volume is assumed: where 2( ) denotes the second invariant of the deviatoric stress tensor = − 1∕3tr( ) and 1( ) = tr( ) the first invariant of the stress tensor. The friction angle is introduced as a function of the equivalent plastic strain and its rate , which will be derived below. This function can be chosen to model the frictional behavior of the shear zone (e.g., based on ring shear test results) and is normally defined through some function ( ) as In order to define the rate of the equivalent plastic strain, the so called spatially rotated rate of plastic deformation is introduced (De Souza Neto et al., 2008): where is the rate of plastic deformation, =̇( ) −1 the plastic velocity gradient and the elastic rotation tensor given by the polar decomposition of the elastic deformation gradient. For the proposed model, the rate of the equivalent plastic strain is defined as A non-associated flow rule is introduced by the Lie derivative of the elastic left Cauchy-Green strain tensor (Simo & Miehe, 1992): with denoting the plastic multiplier. Assuming flow at constant volume, the plastic potential is introduced as By assuming plastic isotropy, it can be shown that the spatially rotated rate of plastic deformation can be written as (De Souza Neto et al., 2008) This type of evolution law emphasizes the analogous extension of the small-strain constitutive law. Inserting Equation 11 into the definition of the rate of equivalent plastic strain (Equation 8), the following evolution law follows:̇=√ 3 In contrast to the over-stress viscoplastic models, the Kuhn-Trucker conditions remain valid and hence the consistency condition completes the definition of the model (W. M. Wang et al., 1997): The model presented here is implemented as a fully implicit exponential return mapping algorithm (Simo & Meschke, 1993). Implementation details can be found in the Supporting Information.

Benchmark Analysis
In a first step, the presented procedure for evaluating co-seismic landslide displacements is applied to the infinite slope conditions shown in Figure 3a for a rate-independent material behavior. The same problem is solved using FEM within the ABAQUS computing environment (Dassault Système Simulia, 2019), providing a benchmark for the MPM simulation. Furthermore, both methods are compared to the Newmark's sliding rigid block analysis to show the influence of compliance within the sliding mass. Infinite slope conditions are incorporated by using periodic boundary conditions. Wave-absorbing infinite elements (Lysmer & Kuhlemeyer, 1969) are added to the bottom of the MPM and FEM model. The only difference between the two procedures lies in the modeling technique of the shear zone. Whereas MPM can deal with the large deformations encountered in the shear zone, special techniques (e.g., remeshing or Eulerian-FEM) would be required for FEM. Therefore, the shear zone is modeled using a contact interface in FEM with the corresponding friction coefficient = tan . The applied parameters for the benchmark analysis are listed in Table 3. The benchmark slope was subjected to the input motion recorded from the Irpinia event in 1980 (RSN 292, H2 direction). The time history was retrieved from the PEER strong motion database (Ancheta et al., 2014).
The evolution of displacements shows that the presented approach of continuous modeling of the shear zone leads to nearly the same results as the corresponding simulation using an interface in FEM (Figure 3b). The slight difference in the displacements is mainly attributed to differences in the numerical scheme (e.g., different shape functions and additional transfer in MPM). Therefore, the FEM and MPM approaches can be seen as identical for this special case of infinite slope conditions, as both account for compliance of the landslide. In contrast, the Newmark's sliding block analysis neglects any wave propagation, but nevertheless it results in rather similar 10.1029/2022JF006618 9 of 25 displacements for the selected parameters. To highlight the influence of compliance, results for MPM are also presented for a softer and stiffer landslide which is reflected in different shear wave velocities v s (Figure 3b). It can be seen that for a lower shear wave velocity the displacements get slightly larger but still remain similar to the other cases. A more detailed discussion on the influence of compliance can be found elsewhere (Kramer & Smith, 1997;Rathje & Bray, 2000). However, the coupled approach proposed by Rathje and Bray (2000) should only be applied with reservation as, owing to the use of a rigid base, seismic waves are trapped in the model, which for most cases results in excessively strong amplifications and an overprediction of displacements.

Model Description
The geometry of the MPM model is based on the cross section presented in Figure 1b, where the shear zone is pre-defined according to the inclinometer measurements with a numerical thickness of = 1.0 m . The groundwater table is introduced at various depths parallel to the surface, and the corresponding water pressure distribution is pre-calculated by solving numerically the problem of saturated flow through porous media. The bedrock is modeled as a linear-elastic base, where the parameters are chosen at the lower end of the range presented in Section 2 as this range covers the bedrock properties up to several hundreds of meters in depth. For the landslide mass, an elasto-plastic model with a Mohr-Coulomb failure criterion is applied using the values for the friction and wave velocities as an average of the presented ranges (Table 4). As shown in the benchmark analysis (Section 3.3), the resulting displacements do not significantly depend on the wave velocities for this range. The use of an average value for the friction angle inside the landslide mass is justified because the range is rather narrow and most of the shearing is concentrated in the shear zone anyway. For the shear zone, the constitutive model from Section 3.2 is applied using different sets of parameters described in the following sections. Owing to the assumed residual state of shearing at constant volume in the shear zone (Section 2.2), the effective pressure of the corresponding material points is kept constant during the seismic analysis and is given by the pre-seismic state.
In the following analyses, the landslide displacements are calculated as the relative movement between a material point from the sliding layer and the base close to the shear zone. The ground motion is applied as both the vertical and horizontal components, where the latter is calculated by rotating the horizontal signals in the direction of sliding.

Effects of Geometry
In the first step, the landslide is investigated without considering any rate-dependent behavior of the shear zone, in order to show the importance of geometrical effects. A strength reduction analysis was performed for an assumed depth of the groundwater of = 4 m to determine the critical friction angle of the shear zone such that equilibrium could just be met. The following seismic analysis is done for different friction angles described by the static safety factor = tan ∕tan using the recorded Norcia input motion.
For an active landslide, the most interesting case is represented for a static safety factor ≅ 1.0 where a Newmark sliding block analysis would suggest that the final co-seismic displacements would tend to infinity. However, the MPM analysis shows, for the Norcia event, a maximal displacement of only 58 mm at the landslide's toe and around 36 mm in the main part (Figure 4a). This illustrates the effect of geometry, where the slide does not move as a single block but is characterized by three zones with different behavior: (a) an unstable upper part, where the friction angle of the shear zone is smaller than its inclination, which represents the driving force of the slide; (b) a stable middle part, which prevents the upper part from a catastrophic failure and decelerates it back to full rest; (c) the toe of the slide, which is not restrained by any soil mass in front of it and therefore is displaced more than the upper part. The latter is particularly amplified by the release of elastic strain energy  stored in this compression zone. Considering the modeled displacements as small to medium, the term geometrical hardening is deliberately avoided.
The comparison to the measurements in borehole A11 shows that the numerical simulation significantly overestimates the landslide displacements. To gain a better understanding of the influence of the strength of the shear zone, the same analysis is performed for different static safety factors (Figures 4b and 4c). This clearly shows that a higher static safety factor leads to smaller landslide displacements, with the corresponding Newmark sliding block analyses (for an average inclination of 8° degrees) showing the fastest decrease. However, even for Newmark, in order to reproduce the measured displacements a safety factor larger than 1.2 would be necessary. Noting that the landslide was slowly moving in the days before the event, such a high safety factor seems unreasonable.

Effects of Rate Dependency
The overestimation of landslide displacement for rate-independent shear zone behavior indicates that it is unreasonable to neglect viscous effects. Moreover, the reported large fluctuation of the groundwater level (Ruggeri et al., 2020) and the creeping velocities emphasize that viscosity should be taken into account. Therefore the critical friction angle of the shear zone was also determined for the minimal (2 m) and maximal (6 m) groundwater level. Based on the continuous monitoring and the direct connection between the landslide movements and the hydrological conditions, a corresponding measured value of velocity was defined for each depth of the groundwater level ( Figure 5). In dry periods, almost no movements were observed, and hence it was assumed that the landslide is exactly in static equilibrium (zero velocity) for a groundwater depth of 6 m with a reference friction angle 0 . To ensure a quasi-static steady state (safety factor = 1 ) for the different water levels, a logarithmic law describing the rate-dependent shear behavior (Equation 6) is introduced as where 0 , and 0 are material parameters. This kind of relation is often chosen for clayey shear zones in active landslides Puzrin & Schmid, 2011;Wedage et al., 1998), but the reference rate 0 is also added in the numerator to avoid the singularity in the logarithm. The application of a logarithmic law ensures that the resistance is not overestimated for higher velocities during a seismic event. The landslide velocity can be linked to the rate of the equivalent plastic strain for a given shear zone thickness , assuming simple shear conditions, aṡ= The same relation can also be applied to link the reference landslide velocity 0 to the reference strain rate 0 . This link between velocity and strain rate (and the same holds for displacement and strain) is critical for defining the thickness of the numerical shear zone, which does not match the real conditions. This procedure is based on the smeared crack approach proposed by Rots et al. (1985) and ensures that the results do not depend on the shear zone thickness. The influence of the viscosity parameter on the friction angle is illustrated in Figure 5 for a wide range of velocities, from the annual creeping rate to the behavior during strong motions. The value of = 0.04 leads to a good fit and will be applied for the following simulations. As has been shown using ring shear testing, a relatively large increase of 4% in resistance per log cycle indicates the presence of a sizable content of clay minerals in the shear zone material (Duong et al., 2018;Scaringi et al., 2018;Tika et al., 1996). For accurate calibration, experimental testing of material from the shear zone for the full range of velocities is necessary. Nonetheless, this serves as a reasonable first assumption for the following calculations and will be investigated in a subsequent sensitivity analysis.
The seismic analysis is carried out for different groundwater levels because the exact depth at the time of seismic events is unknown. This also allows the influence of the pre-seismic conditions on the co-seismic behavior to be investigated. The results illustrated in Figure 6 clearly show that displacements can be modeled more accurately by taking into account rate effects in the shear zone. While for the Norcia event the prediction based on the medium groundwater depth of 4 m provided an excellent fit to the observed permanent deformations, for the Accumoli event the best fit is obtained by assuming = 6 m , and for Visso by assuming = 3 m . This is consistent with the pre-seismic landslide conditions given by the creeping velocity and the precipitation (Figure 1c). The Visso earthquake took place after a rather wet period characterized by increased creeping velocities of around 0.5 mm/day, indicating a higher groundwater level. On the other hand, prior to the Accumoli event, the landslide was moving at a very slower rate, suggesting a lower level of the groundwater.
Another important insight can be obtained by looking at the displacements (1.3 mm) for a groundwater depth of 2 m during the Norcia event. The seasonal fluctuation in the creeping velocity and the measurements of the groundwater level (Ruggeri et al., 2020) suggest that this is what can be expected when a similar earthquake takes place during a period of increased landslide movements.

Assessment of Possible Site Effects
In the previous analyses, the ground motions recorded at the IT-VAL station have been directly applied to the landslide model, as the station is nearby and has similar site characteristics. Ideally, the input motion at the bedrock should be back-calculated using deconvolution from the recorded surface motion for the geophysical profile (defined by e.g., cross-or downhole tests) at the seismic station (Mejia & Dawson, 2006). Such data are not available for the IT-VAL station. However, the Accumoli earthquake was also recorded by two additional seismic stations in the vicinity of the landslide: IT-ASS and IV-ATCC (Russo et al., 2022). The locations of the seismic stations are shown in Figure 1a, and the horizontal ground accelerations (in the direction of the landslide Figure 5. Illustration of the logarithmic rate-dependent friction law for different parameters A (top) and for the calibrated model including the observed data. For groundwater levels d w of 2 m, 4 and 6 m, corresponding critical friction angles, ensuring a quasi-static equilibrium, of 9.5°, 10.8° and 11.8° were calculated by using a strength reduction analysis. Figure 6. Ground motion and evolution of the landslide displacements at AI11 from the material point method analysis for the recorded earthquakes. The displacements are calculated by the difference of the position from a material point above and below the shear zone. The oscillations arise from the elastic waves, which can only be seen because of the small displacements. movement) as well as their Fourier spectra in Figure 7a. The IT-ASS station is located on ground type A (rock) according to Eurocode 8 (Comité Européen de Normalisation, 2004), and the shear wave velocity is estimated as v s,30 = 1,070 m/s (Russo et al., 2022). This makes the IT-ASS record suitable as the input at the bedrock of the La Sorbella model. However, it can be observed that both the horizontal ground accelerations and the Fourier spectra are remarkably similar for the IT-VAL and IT-ASS stations, indicating a lack of significant amplification in the vicinity of La Sorbella landslide. Slightly higher accelerations at the IT-VAL station do not necessarily indicate an amplification: they could also be attributed to the different epicentral distance and topography. In contrast, the  Table 4 and the measured displacement in borehole AI11.
IV-ATCC is located on a 10 m soft layer on top of a stiffer material and is, therefore, classified as ground type E (soft soil) according to Eurocode 8 (Comité Européen de Normalisation, 2004), which is clearly reflected in the signal and its Fourier spectrum (Figure 7a).
The comparison of landslide displacements for the MPM model subjected to the input signals from all the three stations is presented in Figure 7b. Although the final displacements are clearly different, the order of magnitude stays the same and confirms the previous results. One could even conclude that the application of the rock outcrop signal recorded at station IT-ASS leads to a more accurate result, when compared to the measured displacement. The uncertainty in groundwater level, however, makes this conclusion not justified.

Predicted Displacements
The preceding analyses using the ground motions recorded in 2016 allow quantification of the relative effects of geometry and rate dependency on La Sorbella landslide behavior during moderate earthquakes. They also serve as a calibration and demonstration of the presented methodology. Based on this, the behavior for other and especially stronger earthquakes can now be investigated by subjecting the landslide model to a set of recorded ground motions presented in Table 5 (Ancheta et al., 2014), using the same soil parameters as in the previous section. Note that this does not represent a risk assessment for this specific landslide, which would require choosing a representative set of input motions for a defined return period. The goal of this analysis is to get an idea about the order of magnitude for the corresponding landslide displacements, as well as the main factors influencing them.
The evolution of the displacements and the final displacements are shown in Figure 8. As expected, stronger input motions lead to considerably larger displacements which could cause severe damage to local infrastructure, such as the highway in the case of La Sorbella landslide. However, the landslide decelerates quickly and returns to the state of slow motion, as could be expected from the geometrical effects and the rate-hardening constitutive model. Furthermore, since these displacements can be considered as medium to large, geometrical hardening can also contribute to the deceleration, because some weight is transferred from the unstable to the stable part of the slope.
The largest displacements occurred at the toe of the landslide again for all three earthquakes, which can be expected from the unrestrained front due to release of elastic strain energy in the compression zone. However, for the Kocaeli motion the displacement field looks more uniform, which could be due to the impulse-like characteristics of the input signal (Figure 8b).

Influence of Rate Dependency of Shear Strength
So far, the parameters of the shear zone have been calibrated based on the annual creeping behavior of the landslide. This set of parameters also leads to reasonable co-seismic displacements for the recorded input motions. However, these parameters must be applied with caution to stronger earthquakes since the associated displacement rates might be several orders of magnitudes higher. The corresponding constitutive relationship should be calibrated for this range of velocities using appropriate experimental tests instead. The influence of this  Table 5 Earthquake Data for Strong Motions Figure 8. Results from the material point method simulation of La Sorbella landslide subjected to three strong input motions. The evolution of the displacements is shown for different locations alongside the ground motion. For each earthquake, the final co-seismic displacements are presented directly underneath.
rate-dependent behavior is investigated by performing the seismic analysis for different viscosity parameters A, using the Irpinia earthquake motion (Table 5). To ensure quasi-static equilibrium in the pre-seismic state for low viscosities, the reference friction angle is set slightly higher than the critical friction angle ( 0 = 11.0 • ; = 10.8 • ).
In Figure 9a, the evolution of the landslide displacements is presented for both rate-hardening ( 0 ) and rate-softening (  Although a higher viscosity significantly reduces the displacements compared with the rate-independent case ( = 0 ), for the strong Irpinia motion they remain within the same order of magnitude for a reasonable choice of parameter . In contrast, the simulations for the recorded input motions (which were significantly weaker, see Sections 3.3 and 3.4) show a reduction of two orders of magnitude when considering viscosity. This highlights the influence of geometrical effects, which reduce the role of rate effects and are more pronounced for larger displacements. In the case of rate-softening, these effects act against each other, and geometrical hardening prevents the landslide from accelerating into catastrophic failure. The limit case without any geometrical hardening is given by the Newmark's sliding block analysis for an average inclination ( = 8 • ) and the corresponding material model with the same factor of safety ( 0 = 8.15 • ; = = 10.8 • ). The comparison of the final displacement for different viscosity parameters clearly shows that in the absence of geometrical hardening the influence of the rate-dependent behavior is much larger (Figure 9b).

Influence of Potential Softening in the Soil Mass
Little attention has been paid so far to the material of the landslide body and the surrounding soil masses at the bottom and top of the slope. For older active landslides, it is often assumed that the behavior is mainly controlled by the sliding surface where most of the shearing is concentrated. However, a potential strength reduction in the soil mass could counteract the geometrical hardening and might lead to larger landslide deformations. This possibility and its consequences are investigated for La Sorbella landslide, assuming, for stronger motions, linear isotropic softening of the soil mass, where the initial friction angle ( = 30 • ) can reduce to a residual value ( = 16 • ) based on laboratory tests (Ruggeri et al., 2020). The seismic simulation is performed using the Irpinia input motion for two scenarios: (a) Calibrated rate-hardening behavior of the shear zone analogous to Section 3.5 ( = 0.04 ) and (b) rate-softening behavior of the shear zone according to Section 3.6 ( = −0.01 ). The comparison of the displacements for both scenarios and the corresponding results for a constant strength are presented in Figure 10. The consideration of a softening soil mass obviously leads to larger deformations and especially affects the landslide toe. Although the displacements in scenario (a) are significantly higher than for the non-softening analysis, they remain within a certain limit. However, when a rate-softening shear zone is combined with a softening in the soil mass, geometrical hardening does not counterbalance the loss of strength in the shear zone due to the co-seismic acceleration, and the landslide catastrophically fails (Figure 10b).

Geometry and Kinematics
The case study of La Sorbella landslide suggests that the effects of geometry cannot be neglected for calculating accurate co-seismic landslide displacements. This is shown in particular by the rate-independent analysis for different pre-seismic safety factors where, in contrast to the Newmark sliding block analysis, displacements calculated using real geometry do not tend to infinity when the landslide is exactly at equilibrium. However, the extent of these effects is likely to depend strongly on the geometry. La Sorbella seems to be particularly influenced, as the steeper top part is stabilized by a more massive bottom part. This becomes even more evident during strong motions for a rate-softening material in the shear zone, where the transfer of mass from the unstable to the stable part prevents the landslide from accelerating into a catastrophic failure. This effect is likely to be less pronounced in a long landslide that has a more constant inclination of the shear zone.
Based on the monitoring results, Ruggeri et al. (2020) concluded that displacement pattern of the landslide shows only small rotation and very small internal deformations. The inclinometer measurements show a tendency for the upper part to move slightly faster than the landslide's toe. This behavior is observed in constrained landslides to a much greater extent (Oberender & Puzrin, 2016;Puzrin & Schmid, 2012) and is most likely due to the stabilization of the unstable top portion by the stable bottom part. The associated internal deformations are distributed over the entire landslide in the form of elastic deformations leading to increased horizontal stresses in the landslide (Schwager & Puzrin, 2014). In the simulated co-seismic displacements, only small diffusive internal deformation can be observed for large portions of the landslide, even for the stronger motions (Figure 8), which is in agreement with the non-seismic behavior. However, the bottom and top boundaries deform more and for some cases even show distinct internal shearing. The deformation pattern at the top indicates that the location of the shear zone in reality might be different and less localized, with a more complex scarp than the assumed one. In all the simulations, the displacements were largest at the toe, which can be explained by the lack of constraint and the seismically induced release of elastic energy stored in this compression zone. Whereas simulations without considering any rate effects (Figure 4b) show distinct deformations at the toe, more realistic simulations including a rate-hardening shear zone result in a more diffuse pattern, similar to the non-seismic behavior.

Shear Zone Behavior
Based on the yearly displacements and old history of the slide, the total accumulated displacements are likely to be from tens of centimeters to meters. In combination with the sliding in a distinct shear zone, the shear zone is assumed to have not only reached critical state but also the residual shear strength typical for clayey material, undergoing large deformations (Skempton, 1985). Therefore, shearing is expected to occur at constant volume, and the effective pressure can be kept constant during the seismic analysis. Impulsive shearing at fast rates could lead to an increase in void ratio due to a turbulent shear mode (Tika et al., 1996) and lead to a change in pore water pressure. Such behavior would presumably change the post-seismic behavior due to the low permeability as the shear zone contracted, and excess pore pressures would build up. This behavior would be then associated with increased post-seismic displacement rates, which have not been observed for this landslide.
The analysis of La Sorbella landslide also shows that neglecting viscous effects leads to overestimation of co-seismic displacements by two orders of magnitude. Therefore, a logarithmic rate-dependent friction law, calibrated based on the yearly creeping behavior, had to be introduced in the shear zone. This allows an accurate simulation of co-seismic displacements for all three earthquake events recorded in 2016. The presented model can also be used for a risk assessment by subjecting the landslide to other and especially stronger input motions. Therefore, a representative set of earthquake signals for the specific site should be chosen, and the rate-dependent behavior of the shear zone should be calibrated using appropriate experimental testing (e.g., in a ring shear apparatus) for the expected range of velocities. Since proper experimental data are often unavailable, the risk assessment can be carried out in the form of a parameter study to get an idea of the sensitivity and the maximal expected displacement. Whether a proper consideration of rate effects is important for strong input motions depends on how much the landslide is affected by geometrical hardening.

Assessment of Potential Scenarios and Catastrophic Failure
The results presented here suggest that it is unlikely for the active La Sorbella landslide, which is strongly influenced by geometrical hardening, to fail catastrophically. Even a rate-weakening behavior of the shear zone leads to relatively small displacements, and the landslide decelerates back toward a stable state. As this obviously represents a slightly more stable state due to the mass transfer, one might therefore assume that such co-seismic displacements have not happened in the history of the landslide and are less likely to happen again multiple times. However, the upper part of the slide will remain unstable even after several meters of displacement and will still be prone to earthquake-induced displacements. It should be kept in mind that even stable slopes can experience co-seismic displacements (Figure 4).
One of the limitations of the rate-weakening analysis presented here is that it requires the landslide to be in a stable state before the earthquake. Otherwise, the landslide would accelerate without any external influence and would especially not remain in a quasi-static state for a groundwater level varying by 4 m. Therefore, it is not possible to draw any conclusions regarding the landslide's susceptibility to rainfall-induced displacements. In order to investigate this topic, a more sophisticated type of rate-dependent behavior needs to be used to maintain a quasi-static equilibrium and only exhibit softening for higher velocities.
Depending on the pre-seismic conditions and the landslide history, the soil mass downslope from the sliding layer can also be affected by some type of softening such as excess pore pressures and strain localization . The presented analysis based on the example of La Sorbella shows that the co-seismic behavior might be influenced by such a loss in strength and could even fail catastrophically. Considering that La Sorbella landslide was already documented more than 30 years ago and was probably triggered prehistorically in a different topographic situation (Ruggeri et al., 2020), the landslide has most likely experienced large deformations and intensive internal shearing where a part of the potential softening has already occurred. Therefore, catastrophic failure seems rather unlikely, and even for a pessimistic assessment the expected strength loss might be less than in the simulation. The only scenario that leads to a catastrophic failure requires the combination of softening in the landslide mass and in the shear zone. To make a definitive statement on whether such a scenario is plausible, a detailed investigation of the in-situ conditions and experimental testing on undisturbed samples is required.

Application of the Proposed Methodology to Other Landslides
Ideally, the proposed procedure should be tested on further case studies to assess its reliability. Unfortunately, the co-seismic behavior of creeping landslides is sparsely reported, and the required data are mostly unavailable (Bontemps et al., 2020;Lacroix et al., 2014Lacroix et al., , 2015, but we hope that more studies will focus on this important aspect in the future. Accurate simulation of co-seismic landslide displacements also requires high-quality geological and geotechnical models (e.g., geometry, lithological units, constitutive models and the necessary parameters). The proposed methodology is currently limited to plane strain problems, but can be extended to a three-dimensional framework in order to tackle problems with more complex geometry, which, however, will be computationally very expensive.
In order to apply this procedure as a risk assessment to other landslides, the geometry should be reasonably known but can be more complex and incorporate multiple layers and shear zones, since any desired constitutive model can be assigned to each material point. Additionally, an estimate of the geophysical and geotechnical properties of the lithological units is necessary. Based on that, a risk assessment can be made by looking into different scenarios such as different input motions, material models and the corresponding parameters. This will allow identification of the sensitivity to different factors, a probable range of displacements and potential scenarios of catastrophic failure. In cases where experimental data from the shear zone are available, more accurate and sophisticated material models (e.g., Handwerger et al., 2016;Marone, 1998;Ruina, 1983;Scholz, 1998) can be applied, allowing for a more reliable risk assessment. Pore pressure changes controlled by material dilatancy (Iverson, 2005;Iverson et al., 2000), excess pore pressures due to cyclic loading and effects due to frictional heating Vardoulakis, 2000) can also be incorporated. Accounting for these effects would require a hydro-thermo-mechanically coupled MPM, which has been recently applied to a non-seismic landslide simulation by Alvarado et al. (2019), and can be extended to incorporate rigorous seismic boundary conditions. If the coupled approach appears to be computationally too expensive, these effects can be included through the use of phenomenological constitutive models, which account in a simplified manner for the complex underlying physical mechanisms.

Conclusions
An attempt has been made to understand the mechanisms behind the extremely small co-seismic displacements of La Sorbella landslide in Italy during recorded moderate earthquake events, and to investigate implications of these mechanisms for the future landslide behavior during potential strong motions. To achieve this goal, an MPM model is proposed, which includes an accurate representation of the seismic boundaries, ensuring that outgoing waves are not reflected back into the model. Because of the importance of the shear zone and the associated large deformations, a generic large-strain frictional constitutive model for rate-dependent behavior is introduced. The application of MPM brings the advantage that the analysis can be performed regardless of the deformation magnitude. Therefore, the same model can be applied for both moderate earthquakes and stronger motions, which could cause softening and large displacements.
Applying this model to La Sorbella landslide provides several insights: 1. While effects of the rate hardening and of lower groundwater levels on reducing co-seismic displacements are well appreciated, they are not sufficient to explain how the landslide, which was moving up to 0.5 mm/day during certain periods of the year, could experience only 0.4 mm of co-seismic displacement. 2. The answer lies in the geometric effects: weaker ground motions have difficulty in displacing the massive stable middle portion of the sliding body. Note that, owing to the small magnitude of displacements, these favorable geometric effects have nothing to do with the so-called geometric hardening, where soil weight is transferred from the unstable portion into the stable one, which is often the case for stronger motions. 3. Interestingly, the toe of the landslide, which is supposed to be the most stable part because of the upward inclined sliding surface, experiences the largest displacements. This, we believe, is due to the unidirectional unrestraint and the release of elastic energy stored in this compression zone. 4. Rate dependency of strength plays a larger role during light-to-moderate earthquakes. Calibration of the rate-dependent constitutive relationship requires lab tests. However, it has been shown that for weaker ground motions this relationship can also be successfully calibrated based on the observed velocity fluctuations for various groundwater levels during the creeping stage. 5. For stronger motions, co-seismic displacements become larger, causing a weight transfer between unstable and stable parts, and activating geometric hardening. As a result, even the rate-softening of strength sometimes observed in the shear zone would not bring the landslide to catastrophic failure. Such a failure would become possible only if the soil mass in the sliding body and further downhill experienced strain softening, which is unlikely in the case of La Sorbella landslide.
The study presented here demonstrates that the co-seismic behavior of La Sorbella landslide is strongly influenced by the interplay of geometrical and viscous effects, as well as potential softening in the shear zone and the surrounding soil mass. As a next step, it would be useful to investigate effects of some typical slide geometries on the seismic behavior of slow-moving landslides. As a possible broader impact, we hope that the proposed approach could pave a way toward more reliable risk assessment for slow-moving landslides.

Data Availability Statement
Data regarding La Sorbella landslide have been taken from Ferretti et al. (2019) and Ruggeri et al. (2020). Precipitation data from the Casanuova dam weather station are accessible at the web portal of Servizio Idrografico Regione Umbria (https://annali.regione.umbria.it). All data regarding the seismic stations and the recorded earthquake motions (Accumoli, Visso and Norcia) processed by Russo et al. (2022) are accessible at the ITACA web portal (https://doi.org/10.13127/itaca.3.2). The input signal for the strong motions has been retrieved from the PEER Ground Motion Database NGA-West2 (Ancheta et al., 2014) and is accessible at the PEER web portal (https://ngawest2.berkeley.edu/). The digital terrain models are provided by Ministero dell'Ambiente e della Tutela del Territorio e del Mare and are accessible at their web portal (http://www.pcn.minambiente.it/).