Numerical analysis of downward progressive landslides in long natural slopes with sensitive clay

Landslides occurring in sensitive clay often result in widespread destruction, posing a signi ﬁ cant risk to human lives and property due to the substantial decrease in undrained shear strength during deformation. Assessing the consequences of these landslides is challenging and necessitates robust numerical methods to comprehensively investigate their failure mechanisms. While studies have extensively explored upward progressive landslides in sensitive clays, understanding downward progressive cases remains limited. In this study, we utilised the nodal integration-based particle ﬁ nite element method (N-PFEM) with a nonlinear strain-softening model to analyse downward progressive landslides in sensitive clay on elongated slopes, induced by surcharge loads near the crest. We focused on elucidating the underlying failure mechanisms and evaluating the effects of different soil parameters and strain-softening characteristics. The simulation results revealed the typical pattern for downward landslides, which typically start with a localised failure in proximity to the surcharge loads, followed by a combination of different types of failure mechanisms, including single ﬂ ow slides, translational progressive landslides, progressive ﬂ ow slides, and spread failures. Additionally, inclined shear bands occur within spread failures, often adopting distinctive ploughing patterns characterised by triangular shapes. The sensitive clay thickness at the base, the clay strength gradient, the sensitivity, and the softening rate signi ﬁ cantly in ﬂ uence the failure mechanisms and the extent of diffused displacement. Remarkably, some of these effects mirror those observed in upward progressive landslides, underscoring the inter-connectedness of these phenomena. This study contributes valuable insights into the complex dynamics of sensitive clay landslides, shedding light on the intricate interplay of factors governing their behaviour and progression.


Introduction
Landslides in regions characterised by sensitive clay are a recurring concern, particularly prevalent in northern regions of the world.These events pose severe threats to human life and can cause significant economic losses due to their widespread and devastating impact.As strain accumulates, the shear strength of sensitive clay decreases (Pusch, 1966), a phenomenon known as strain-softening.Sensitivity, defined as the ratio of the peak to residual shear strength, serves as a fundamental parameter that signifies the degree of strain-softening that clay can undergo and plays a pivotal role in the failure mechanisms of geostructure within sensitive clay (Skemption and Northey, 1952).For soils exhibiting high sensitivity levels, they lose structural stability and undergo a transformation from a solid state to a liquid-like substance due to the significant decrease in shear strength as strain increases.This phenomenon is termed the remoulding of sensitive clays, and the residual shear strength associated with this process is referred to as remoulded shear strength.In a typical retrogressive landslide occurring in sensitive clays, the clay located adjacent to the initial shear band within a slope shifts from its original position due to strain-softening behaviour.This subsequently gives rise to instability within the newly formed back scrap.Consequently, a succession of failures, commonly referred to as progressive failures, i.e. unfolds, which ultimately lead to widespread destruction (Bjerrum, 1955).Despite the numerous numerical methods developed for modelling landslide problems involving large deformation, there has been limited research dedicated to analysing postfailure scenarios of landslides in sensitive clay (Urmi et al., 2023).This underscores the pressing necessity to both employ and advance novel numerical tools to attain a more profound comprehension and predictive capacity pertaining to these catastrophic events.
Considering the failure mechanisms of sensitive clay landslides, three primary types stand out: flow slides (or earth flows) (Fig. 1a), translational progressive landslides (Fig. 1b), and spreads (Fig. 1c) (Locat et al., 2011).Additionally, in situations of lower severity, such as in cases of low sensitivity or slow strain-softening rate, only a single flow slide is typically observed.Progressive failures can be classified as upward and downward failures (Dey et al., 2016a), or uphill and downhill failures (Bernander et al., 2016), depending on the location of the disturbing agents.Upward progressive failure is commonly referred to as 'retrogressive failure', while downward progressive failure is termed 'progressive failure' in certain studies (Quinn et al., 2012).In upward progressive failure scenarios, shear bands or local failures typically begin near the slope's toe and extend towards areas further upslope.In contrast, in cases of downward progressive failure, shear bands may start near the loaded area positioned upslope and then spread downward, resulting in progressive failures.
Flow slides in sensitive clays, initiated by erosion over time, start as single rotational failures, often progressing upward, as noted by Bjerrum (1955).Translational progressive landslides, moving as whole flake slides along elongated surfaces parallel to the ground (Varnes, 1958), are common on lengthy, gently sloping terrain, especially in a downward direction, typically triggered by surcharge loads near the slope crest (Bernander et al., 2016;Tremblay-Auger et al., 2022;Zhang et al., 2023a).Combined failure of flow slides and translational progressive landslides may occur due to excavation at the foot of a long natural slope, as observed in the progressive landslides in the Xiangjiaba reservoir area, Southwest China (Xu et al., 2014).In spreads, weak soil layers undergo stress-induced horizontal failure, resulting in soil mass displacement and the formation of triangular soil blocks called horsts and grabens (Mollard and Hughes, 1973).
Later, several numerical approaches were devised and adopted to examine the post-failure dynamics of landslides occurring in sensitive clay.These methods include the coupled Eulerian-Lagrangian (CEL) approach (Islam et al., 2019;Wang et al., 2022), the material point method (MPM) (Wang et al., 2016;Tran and Sołowski, 2019), the particle finite element method (PFEM) (Zhang et al., 2017(Zhang et al., , 2020)), the smoothed particle finite element method (SPFEM) (Yuan et al., 2020;Zhang et al., 2018), and the remeshing and interpolation technique with small strain (Shan et al., 2021), among others.Most of these investigations have centred on upward landslides, particularly those associated with flow slide failure or spread failure resulting from toe erosion on the slope (Locat et al., 2013;Yuan et al., 2020;Zhang et al., 2020;Wang et al., 2022), or seismic loading (Islam et al., 2019).When surcharge loads and erosion are considered together (Dey et al., 2016a, b), a combined upward and downward progressive failure occurs.This can also be observed in sensitive clay on a steep and long slope subjected solely to gravity loads (Wang et al., 2016).Despite these contributions, there remains a scarcity of research dedicated to the comprehensive analysis of purely downward landslides occurring on extensive natural slopes composed of sensitive clay.The processes governing the failure of downward progressive landslides after they occur, as well as the influence of material characteristics on shaping this phase, remain areas where our understanding is incomplete.Given these considerations, the current investigation seeks to explore the failure mechanism within sensitive clay on a long slope affected by surcharge loads near its peak, with no consideration of erosion at the slope's base.The strength of sensitive clay diminishes nonlinearly as deviatoric plastic strain increases.Simulations are performed using the nodal integration-based PFEM (N-PFEM) (Meng et al., 2021;Zhang et al., 2022).The investigation focuses on the impact of four pivotal aspects: the thickness of the sensitive clay layer, the strength gradient along the depth of clay, the clay sensitivity, and the rate of strainsoftening.These variables are systematically investigated to determine their effects on failure modes and the extent of diffused areas.
The subsequent sections of the paper are structured as follows.Section 2 outlines the problem under consideration.Section 3 introduces the formula used to describe the strain-softening feature of sensitive clay.The simulation method is summarised in Section 4 with the model setup presented in Section 5. Section 6 shows the simulation results with accompanying discussion.Finally, Section 7 draws conclusions based on the findings.

Problem description
To investigate the downward landslide in sensitive clay, we consider a typical model, as illustrated in Fig. 2. The model consists of two sections: one extending over a horizontal distance of 190 m to the right with an inclination of 5:71 , and the other spanning 100 m to the left with no slope.The total depth of the model, denoted as H, is 20 m.The bottom, left and right boundaries are fixed in simulations.
Due to the long-term weathering processes, a non-sensitive crust layer typically forms at the top, exhibiting higher undrained shear strength compared to the layer below it.The thickness of this crust layer generally varies from 1 m to 6 m (Lefebvre et al., 1987).In this study, we assume a uniform depth, H c , of 5 m for this crust and a constant shear strength of 50 kPa.Below the 5 m crust layer, a clay layer extends to a depth of 15 m.Although the average shear strength of the clay remains steady at 40 kPa, its shear strength escalates proportionally with depth according to the following equation: where s u0 represents the initial shear strength of the clay (kPa), s ug is a constant (kPa), k denotes the strength gradient (kPa/m), and z represents the depth of the soil below the ground surface (m).
For the clay layer with linearly increasing shear strength, there is a sensitive clay layer at the base of the thickness H s , and strainsoftening characteristics are applied to this layer.The unit weight, g, of the crust and clay layers is the same, which is 18 kN/m 3 .To induce failure, an embankment measuring 20 m in width and 7.5 m in height is positioned on the illustrated model.The embankment is situated 150 m away from the slope toe.The unit weight of the embankment is represented by g E .

Strain-softening representative of sensitive clay
For the sensitive clay layer, the evolution of shear strength (cohesion), s u , is governed by established softening rules that correlate it with the deviatoric plastic strain (ε p d ).As ε p d increases, s u progressively decreases until it reaches its residual value.Numerous experimental tests have consistently demonstrated a highly nonlinear relationship between the decreasing shear strength (s u ) and the increasing deviatoric plastic strain (ε p d ).A comprehensive summary of these findings has been provided by Urmi et al. (2023).
In this study, we describe strain-softening using mathematical equations from the work of Yerro et al. (2016): where s up and s ur represent the peak and residual values of the cohesion, respectively.The softening rate is regulated by h.A larger magnitude of h represents a faster strain-softening.The deviatoric plastic strain invariant ε p d is defined as where e p ij is the deviatoric part of the plastic strain tensor.To quantify the degree of strength reduction, sensitivity is adopted (Skemption and Northey, 1952): which has been widely used for characterising sensitive clay (Dey et al., 2015(Dey et al., , 2016b;;Tran and Sołowski, 2019;Zhang et al., 2020).

Simulation methodology
In recent decades, there has been notable progress in developing various numerical methods and using them to tackle large deformation geotechnical problems.Of these methods, PFEM was initially introduced to simulate the flow of free surface evolution and fluidestructure interaction (Idelsohn et al., 2004(Idelsohn et al., , 2006)).It combines particle-based methods and the traditional Lagrangian finite element method, allowing it to handle large deformation problems accurately.Since its invention, the PFEM has proven its efficacy and robustness in simulating various challenging scenarios, including fluidesolid interaction, multiphase flow, and problems involving phase change such as melting.However, when applied to geotechnical problems, the classical PFEM has encountered some challenges.
One common challenge associated with the PFEM is volumetric locking, which occurs when using three-node triangular elements to simulate elastoplastic materials such as clay.To address this issue, high-order elements were implemented in a modified version of the PFEM (Zhang et al., 2013).However, this approach necessitates variable mapping for both quadrature points and mesh nodes when handling history-dependent materials.Another solution is to adopt a mixed finite element formulation to overcome the volumetric locking issues (Monforte et al., 2017(Monforte et al., , 2018)).This allows the use of low-order elements such as three-node triangular elements.However, information still needs to be transformed from the old meshes to the new meshes when simulating history-dependent materials, and errors can accumulate during the variable mapping process.The smoothed PFEM (SPFEM) (Zhang et al., 2018), which involves calculating and storing the information of all field variables at nodes, was proposed to eliminate the requirement of variable mapping issues.Notably, a temporal instability issue arises when the smoothing technique is adopted in the conventional displacement finite element framework.Consequently, ad-hoc stabilisation techniques are needed to prevent stress oscillations in dynamic analysis (Jin et al., 2021;Shafee and Khoshghalb, 2022).
The N-PFEM is adopted in this study, which is an enhanced PFEM.In the N-PFEM, a mixed variational principle, specifically, the generalised Hellinger-Reissner variational principle (Zhang et al., 2019), is implemented within a nodal integration scheme as part of the finite element formulation.Consequently, variable mapping is eliminated, and information is computed and stored in particles.Moreover, the N-PFEM exhibits natural temporal stability when employing an implicit time integration, eliminating the need for ad-hoc stabilisation techniques and enabling the use of larger time steps to enhance computational efficiency.
In the N-PFEM, integration is performed on cells rather than finite elements, as illustrated in Fig. 3. To construct these cells, the centroid and three mid-edge points of each triangle are connected, resulting in three quadrilaterals of equal area (equivalent to onethird of the triangle area).The integrated cell for each particle is formed by aggregating all one-third portions of the adjacent triangles.The N-PFEM does not require variable mapping since all information is stored and calculated directly on particles.The computational cycle of the N-PFEM is described as follows: (1) Use a cluster of particles to define the domain (Fig. 4a); (2) Perform Delaunay triangulation with the particles to establish the boundaries of the domain (Fig. 4b); (3) Construct cells through triangular meshes, as described in Fig. 3, and subsequently solve the equations governing the problem (Fig. 4c); (4) Update particle positions and erase old meshes (Fig. 4d); (5) Repeat Steps (2)e( 4) for all time steps.The adopted N-PFEM has undergone rigorous validation for large deformation problems through simulations across various scenarios.These include granular column collapse under both quasi-static (Zhang et al., 2022) and dynamic conditions (Zhang et al., 2023b), underwater granular flow and induced tsunamis (Zhang et al., 2022), cliff recession (Meng et al., 2021), deep penetration, and plate anchor pullout (Wang et al., 2023), among others tests.Mesh dependency may arise when modelling strainsoftening behaviour, owing to ill-posed boundary-value problems in classical plasticity theories (Bazant et al., 1984).Several regularisation techniques address the mesh dependency of strainsoftening behaviour by introducing a length scale into boundary value problems, either explicitly or implicitly, such as the ones using viscoplastic models (Wang et al., 1997;Zhang et al., 2017;De Borst and Duretz, 2020), gradient plastic models (De Borst andMühlhaus, 1992), andCosserat continuum (De Borst, 1991;Huang and Bauer, 2003;Tang et al., 2018;Wang et al., 2019).A recent study by Zhang et al. (2023b) on the mesh dependency effects in modelling landslides in sensitive clay showed that larger meshes can result in a larger diffused area and may affect collapse times, but not the failure modes.Thus, although simulation results are associated somewhat with mesh dependency, the findings in failure modes are still valid.

Model setup and simulations
This study employs the Tresca yield criterion to conduct rigidplastic dynamic analysis.Specifically, the clay layer subjected to an additional consideration regarding sensitivity, with a specific layer at its base having a thickness denoted as H s , exhibits strainsoftening.The model was initially discretised using meshes sized at 0.6 m (0.03H).This yielded a total of 39,817 meshes and 20,437 particles.The simulation employed a time step of 0.01 s.To counteract the impact force arising from the abrupt placement of the embankment, we introduced a gradual increase in the unit weight of the embankment over the first 10 time steps, spanning a duration of 0.1s.The initial stress condition is set by the surcharge load from the embankment.This research explores four specific factors and their influence on the processes of downward landslides: (1) the thickness of the sensitive clay layer (H s ), (2) the strength gradient (k), (3) the sensitivity (S t ), and ( 4) the strain-softening rate parameter (h).All simulations were executed on a desktop computer outfitted with a 3.7 GHz CPU and 16 GB of memory, operating on Microsoft Windows Server (Version 10.0).The optimisation solver MOSEK (Version 10.1) within the MATLAB environment (R2023a) was employed to solve the resulting problems, as per our previous study (Zhang et al., 2023b), ensuring computational efficiency.For instance, the computational time for Case 4 shown in Table 1 amounted to 12 h for 1740 analysis steps.
The linear increase in undrained shear strength with depth in sensitive clay has been tested in field investigations (Low et al., 2010;Locat et al., 2017).For instance, a strength gradient (k) of 2.66 kPa/m was found in Champlain Sea clays (Wang, 2019), and 2.5 kPa/m in Sainte-Monique (Locat et al., 2015).According to the sensitivity classification by Skemption and Northey (1952), the sensitivity (S t ) ranges from 1 to 2 in low sensitive clays, 2 to 4 in medium sensitive clays, 4 to 8 in highly sensitive clays, greater than 8 in extra-sensitive clays, and greater than 16 in quick-clays.Regarding plastic strain leading to residual strain-softening in sensitive clays, field test results have provided various values, including 0.125e0.52 in Eastern Canadian (Locat et al., 2013), 0.55 in Sainte-Monique, Quebec (Locat et al., 2015), 2 in Saint-Barnab e-Nord, and 1.8 in Saint-Jude, as determined by extrapolation from the direct simple shear (DSS) curve (Urmi et al., 2023).
Based on the observations mentioned above, we obtain the following material parameters to describe the softening feature.The strength gradient (k) in Eq. ( 1) is considered within the range of 1.2e2.8kPa/m.The sensitivity (S t ) is considered to be between 2 and 4, and the strain-softening rate parameter (h) is varied from 5 to 20.Fig. 5 provides an illustration of the evolution of undrained shear strength (s u ) for a case with S t ¼ 2 and s u0 ¼ 40 kPa in the strain-softening model.To achieve 99.9% degradation of shear strength, the corresponding values of deviatoric plastic strain are 1.4,0.7, 0.5, and 0.35, respectively, with the strain-softening rate parameter (h) being 5, 10, 15, and 20.Table 1 presents a summary of cases, including detailed model set-up information and the final collapse results in the post-failure process.It should be noted that Case 1 is simulated in an insensitive scenario, exhibiting a weak single-flow slide, and achieving stability only after 6.9 s.

N-PFEM simulation results
To demonstrate the significant influence of sensitive clay on the extent of diffusion in landslides, we first simulate an insensitive case, Case 1, with all soil parameters detailed in Table 1.In this case, the applied unit weight of the embankment (g E ) is set to its minimum value (27 kN=m 3 ) to induce slope instability.Fig. 7 shows the final deposit at time t ¼ 6.9 s with the deviatoric plastic strain (ε p d ) mapping and the displacement (u) mapping revealing the shear band and the degree of damage, respectively.As the displacement due to surcharge loads from the embankment increases, zones with higher deviatoric plastic strain (ε p d ) indicate the failure surfaces, suggesting a flow slide (Fig. 7a).However, this flow slide leads to limited deformation and stabilises quickly.There are no progressive collapses since the undrained shear strength remains unaffected by the increasing ε p d .Therefore, the failure mode in Case 1 is a single flow slide.The final diffused area is confined primarily near the embankment, with a final diffused distance (L D ) of 44 m (Table 1).It is worth noting that although L D is non-negligible, the damage degree at the ground surface is not substantial, and the final embankment displacement (U E ) is merely 0.5 m (equivalent to only 2.5% of the upper width of the embankment), as shown in Table 1.Furthermore, the displacement within the diffused area is quite limited, with the maximum value reaching just 0.677 m (Fig. 7b).
Next, we simulate a scenario with a very thin sensitive clay layer (H s ¼ 2 m) at the bottom of the slope, which is presented as Case 2 in Table 1.In this case, the minimum value of g E is the same as that in Case 1.However, the scale of collapse is significantly larger than that in Case 1, exhibiting a combination of a single flow slide and a translational progressive landslide.Fig. 8 illustrates the collapse process of Case 2 until stabilisation, using the deviatoric plastic (1) SF e Single flow slide; TR e Translational progressive landslide (e.g.Fig. 1b); FL e Flow slides (e.g.Fig. 1a); SP e Spread (e.g.Fig. 1c).
(    8a), with a shape and location similar to that in Case 1. Subsequently, the width of the failure surface caused by this single flow slide progressively widens, giving rise to two shear bands within it.Furthermore, there is a noticeable uplift on the ground surface at time t ¼ 6.9 s (Fig. 8b), which, in contrast, is the time that the slope stabilises in Case 1.During this period, inclined shear bands emerge between the flow slide and the embankment.Following this, a translational progressive failure, characterised by slow downward extension, becomes evident within the sensitive clay layer at the bottom (Fig. 8c).In this process, the soil block above the translational progressive failure collapses as a whole entity and extends downwards.The translational progressive failure halts at the foot of the slope, resulting in ground surface heaving (Fig. 8d).
As indicated in Table 1, the final diffused distance (L D ) measures 158.8 m in Case 2 (e.g. with a thin sensitive clay layer), which is 260.9% larger than that observed in Case 1 (e.g.no sensitive clay layer).Additionally, the embankment displacement (U E ) reaches 6 m, equal to 30% of the upper width of the embankment.Apparently, the simulated results underscore that even an exceedingly thin sensitive clay layer can trigger a large-scale collapse in a slope.
Cases 1 and 2 illustrate two typical failure modes: a single flow slide, and a progressive landslide combined with a single flow slide and a translational progressive collapse.Through a series of simulations by varying the model setup in this study, it was observed that in a long slope, multiple flow slides or spreads may also occur within a progressive landslide.The mechanism and degree of collapse are influenced by four key factors: the thickness of the sensitive clay layer (H s ), strength gradient (k), sensitivity (S t ), and strain-softening rate (h).The upcoming content will delve into the impact of these four factors, primarily presenting results using mapping figures with the deviatoric plastic strain (ε p d ), displacement (u), or velocity (v).To distinguish between different failure mechanisms, we will use the notations F, T, and S (with corner markers to indicate the collapse order) to represent flow slides, translational progressive landslides, and spreads, respectively.

Effects of sensitive clay thickness
Cases 2e4 in Table 1 present three simulations with sensitive clay thickness (H s ) at the bottom of the slope being 2 m, 6 m, and 15 m, which give the applied minimum unit weight, g E , of embankment required to cause collapse as 27 kN/m 3 , 26 kN/m 3 , and 22 kN/m 3 , respectively.All three cases are simulated with a strength gradient (k) of 1.2 kPa/m, displaying that the initial undrained shear strength (s u0 ) linearly increases with increasing depth, from 31 kPa at the top of the layer below the crust to 49 kPa at the bottom.Fig. 9 shows displacement mapping in the final deposits, and Fig. 10 presents velocity mapping in the first 10 s for Cases 2e4.The peak values represent the maximum displacement or velocity in each figure.When collapses start under surcharge loads, all failures indicate a single flow slide (F 1 ) near the embankment, exhibiting the same failure mode as in Case 1 with insensitive clays (Fig. 10aed, and g).The shape of F 1 is not affected by the sensitive clay thickness (H s ).
After F 1 occurs, with an increasing sensitive clay thickness, the failure mechanism transitions from a translational progressive landslide to progressive flow slides.The remoulded clay flows underneath and squeezes up the interface between the insensitive clay layer and the sensitivity clay layer, forming a long shear band (T 1 ) in a translational progressive landslide (Fig. 10c and f).Meanwhile, several collapses, namely F 2 , F 3 , and F 4 , quickly appear in downward progressive flow slides (Fig. 10h and i), but F 3 and F 4 only arise from the middle section of collapse.In Cases 2 and 3, which exhibit translational progressive failure, there is a noticeable difference in the affected areas (Fig. 9a and b).Due to the thinner sensitive clay, T 1 in Case 2 (Fig. 10c) forms a smaller angle with the slope bottom compared to T 1 in Case 3 (Fig. 10f), resulting in a longer diffused distance, L D , of 158.8 m in Case 2 compared to 92.1 m in Case 3 (Table 1).To exclude the effects of g E , Case 5 is simulated with the same g E value (27 kN =m 3 Þ as in Case 2. Simulation results show that increasing g E from 26 kN/m 3 to 27 kN=m 3 results in a 5.86% increase in L D , indicating a limited influence.Therefore, when comparing Case 2 (H s ¼ 2 m) to Case 5 (H s ¼ 5 m), with all other parameters being the same, it is observed that a thinner sensitive clay layer causes a larger diffused area when the failure indicates a translational progressive landslide.
According to velocity (v) mapping in Fig. 10, the collapse velocity is always slow (v < 2 m/s) in a translational progressive landslide, while it increases rapidly in progressive flow slides, exceeding 5 m/ s within just 1 s (Fig. 10g).However, even though the collapse speed remains relatively low compared to progressive flow slides, a translational progressive landslide could still result in significant damage due to an extremely long collapse time and large-scale areas affected.For example, Case 2 takes 73.6% more time than Case 4 but covers only 13% less diffused distance than Case 4 (see Table 1).In addition, velocity mapping is smoothed without obvious differences between regions in a translational progressive landslide, thus the beginning of collapse (i.e.areas near the embankment) may also take a considerable amount of time to stabilise (e.g.Case 2, as shown in Fig. 10aec).However, areas near the embankment will stabilise earlier in a progressive downward failure of flow slides, as the velocity is primarily concentrated in the middle and end sections of the landslide after F 3 has completely formed (Fig. 10i).Additionally, a short and quickly stopping spread failure may appear from each flow slide, expressed as S 1 , S 2 , and S 3 , respectively, propagating underneath the crust layer (Fig. 10h and  i).

Slope with a thin sensitive clay layer
In the case of very thin sensitive clay (H s ¼ 2 m), Cases 6e9 in Table 1 provide four typical examples of two failure mechanisms: (1) single flow slide (SF), and (2) combined single flow slide (SF) and translational progressive landslide (TR), as demonstrated in Fig. 11, with deviatoric plastic strain (ε p d ) mapping in the final deposits.In the range of 5e20 for the strain-softening rate parameter (h), h does not change the failure mode, and even its effects on diffused areas are quite small.However, an increase in sensitivity (S t ) may lead to the transformation of a single flow slide to a combined single flow slide and translational progressive landslide (see Fig. 11).Additionally, a higher strain-softening rate parameter (h) tends to result in a rapidly formed clear translational progressive shear band (T 1 ) (Fig. 11d).In this scenario, the first flow slide (F 1 ) is relatively vague.Furthermore, it is clear that no obvious straight inclined shear bands (expect F 1 ) emerge over the large area of the landslide (Fig. 11c and d), since the landslide collapses as a whole section within a large-scale plastic zone above the failure surface in a translational progressive failure were reported by Dey et al. (2016b).
Additional simulations are performed to analyse the effects of sensitivity (S t Þ, ranging from 1 to 4, on the final diffused distance (L D Þ and the final embankment displacement (U E ), as indicated in Fig. 12.There is a demarcation line, S t ¼ 1:9, that distinguishes  between different failure mechanisms.If the sensitivity value is lower than the demarcation value, it indicates a single flow slide.Conversely, if the sensitivity value is higher than the demarcation value, it suggests a combined occurrence of a single flow slide and dominant translational progressive landslide.However, this demarcation value (S t ¼ 1:9) may be influenced by several factors, making it stringent and applicable primarily under conditions where k ¼ 1.2, h ¼ 10 and H s ¼ 2 m.Final U E increases almost linearly with an increase in S t within the range of 1e4, while final L D only demonstrates an almost linear relationship with S t in cases dominated by translational progressive landslides (1:9 S t 4Þ.Furthermore, the growth rate of U E is significantly less than that of L D with an increase in sensitivity in a translational progressive failure.Regarding the effects of the strain-softening rate parameter (h), Fig. 13 illustrates the evaluation of L D and U E over time.
Although within the range of 5e20, h has a relatively minor impact on the extent of damage, it does affect the time required for stabilisation.Fig. 13 indicates that more time is required to re-stabilise with a smaller h.For example, collapse in Case 8 (Fig. 11c) takes 49.2% more time to stabilise compared to Case 9 (Fig. 11d), while the difference in diffused distance is only 1.7% less than that in Case 9.In addition, as mentioned earlier, in a translational progressive failure, speeds remain similar from the beginning to the end of the landslide, resulting in the same stabilisation time for all sections of the collapse.

Effects of strength gradient
Based on the previous simulations, a local shear band, indicating a flow slide located near the embankment, is certain to appear first if collapse starts.Its shape and size remain unaffected by the sensitive clay thickness (H s ), sensitivity (S t ), and strain-softening rate parameter (h).In contrast to the effects of these three parameters, the strength gradient (k) significantly affects the size of the first shear band.Fig. 14a presents the sketch and dimensions of the first shear band for different strength gradients (k) in the range of 1.2e2.8kPa/m.The radius of the first shear band decreases with increasing k, and both the span and depth decrease almost linearly with increasing k (Fig. 14b).A comparable phenomenon, wherein an increase in k reduces the size of shear bands, has been noted in upward progressive flow slides (Wang et al., 2022), potentially influencing the post-failure process of all such slides.However, it seems that k does not impact the sizes of the second flow slide (F 2 ) in a downward progressive landslide, which always passes through the bottom area of the slope.Fig. 15 shows the final deposit with deviatoric plastic strain (ε p d ) mapping for Cases 10, 11, and 4 in Table 1 with three different values of k, i.e. 2.4 kPa/m, 1.6 kPa/m, and 1.2 kPa/m, respectively.It is evident that a smaller strength gradient will result in a longer diffused distance in downward progressive landslides, mirroring a similar effect observed in upward progressive landslides (Wang et al., 2022).
Fig. 16 presents velocity (v) mapping in the first 7.5 s of the collapse for Cases 10, 11, and 4. As the strength gradient (k) decreases, the rate of failure increases, potentially leading to a higher number of collapses.For instance, when k equals 2.4 kPa/m, 1.6 kPa/ m, and 1.2 kPa/m, the flow slides occur 2, 3, and 4 times, respectively, with each collapse initiating earlier as k decreases (Fig. 16).Starting from the third flow slide (F 3 ), collapse only commences from the middle section to the end section of the landslide (Fig. 16feh and i), resulting in a notable reduction in speed for the area near the embankment.Simulated outcomes suggest that the slope base serves as a critical pathway for the second flow slide (F 2 ) following F 1 (Fig. 16bed and g).A smaller value of k corresponds to a lower undrained shear strength at the slope base.Therefore, with a reduced strength gradient, the likelihood of occurrence of downward progressive flow slides increases, along with the severity of resultant damage.Although spread failure may appear from each flow slide, spread failure stops very quickly because of the quick progressive flow slides and the limited distance between the adjacent two flow slides.Therefore, some parallel short shear bands appear in a dense distribution in the crust layer (Fig. 15).
The above simulations present a larger moving soil block produced by the first flow slide with the decreasing strength gradient (k), causing a more drastic and quicker collapse.Therefore, when k is high in a case with a very thin sensitive clay layer at the bottom of the slope, indicating a higher initial undrained strength at the bottom of the slope, the model is relatively more stable.For example, further examples indicate that when k is taken as 1.6 kPa/ m in Case 2 (k ¼ 1:2 kPa=m), the final diffused distance is 155.5 m, which is less than 2.1% of the final diffused distance of 158.8 m in Case 2. When k is taken as 2 kPa/m in Case 2, the model is stable.
In cases of low strength gradient combined with a deep sensitive clay layer, a translational progressive landslide may occur if the sensitivity or strain-softening rate is high.Figs. 17 and 18 illustrate two typical examples, one with high sensitivity (Case 12) and the other with a high strain-softening rate (Case 13), both at t ¼ 3.5 s.In these examples, the strength gradient (k) is set to 1.2 kPa/m, and sensitive clays are present beneath the crust layer.In situations  with a low strength gradient, the potential for a dire collapse increases when there is high sensitivity or a rapid strain-softening rate.In such cases, a long shear band, denoted as T 1 , can develop after the initial flow slide (F 1 ).This shear band extends across the slope base and terminates at the slope's foot due to the diminished undrained shear strength at the bottom (Fig. 17a and 18a).The failure caused by F 1 signifies a translational progressive landslide.However, unlike slowly evolving failure involving a thin sensitive clay layer, this type of collapse occurs at a high velocity (Fig. 17b and  18b), encompassing the entire depth of the slope towards its end.This translational progressive failure moves as a single flake slide because the long failure surface at the bottom of the slope forms instantaneously, as reported by Varnes (1958).Additionally, it is worth mentioning that some inclined shear bands may appear (e.g.Fig. 17a) since the landslide does not collapse within a wholly plastic zone, owing to the combined effects of F 1 and T 1 .To conclude, the combination of a low strength gradient and deep sensitive clays with high sensitivity or a high strain-softening rate poses a very dangerous scenario.

. Effects of sensitivity
To focus on the downward progressive failures of sensitive clay landslides, we will use a slope with a high strength gradient (k) set at 2.8 kPa/m to discuss the effects of sensitivity (S t ) and the strainsoftening rate parameter (h).Cases 14e16 in Table 1 are three simulations with sensitivity (S t ) of 1.5, 2, and 4, all considering sensitive clay thickness of H s ¼ 15 m.However, Case 14 is simulated for the first 15 s due to boundary limitations.
Fig. 19 illustrates the deviatoric plastic strain (ε p d ) mapping for these three cases during the initial 15 s.A higher sensitivity value (S t ) intensifies the collapse process, resulting in a greater diffused distance.Moreover, the simulated results reveal that in downward landslides involving sensitive clay, spread failure is more likely to occur with a low sensitivity value (e.g.Fig. 19c), whereas high sensitivity can transform spread failure into progressive flow slides (e.g.Fig. 19i).Such effects mirror the observations in upward landslides with sensitive clay (Dey et al., 2015;Wang et al., 2022).Fig. 20 illustrates the collapse process in Case 16, displaying the displacement (u) mapping during the first 15 s.In this scenario, a total of 5 flow slides crop up during the first 15 s, indicating a significant and serious flow slide failure.Due to the higher shear strength in the crust layer compared to the upper layer of sensitive clays, the crust moves more slowly than the clays beneath it during the failure, leading to a misalignment between these two layers (e.g.Fig. 20a).Notably, there is no evident spread failure expanding from the first flow slide, F 1 , as the second flow slide, F 2 , occurs very rapidly (Fig. 19ge20a).In contrast to the third flow slide, F 3 , in Case 11 (see Fig. 16f) and Case 4 (see Fig. 16h), where the strength gradient (k ¼ 1.6 kPa/m in Case 11 and k ¼ 1.2 kPa/m in Case 4) is lower and it originates solely from the middle section of the landslide, F 3 in Case 16 (k ¼ 2.8 kPa/m) commences near the top of the slope, traverses the slope base, and extends acrosses the entire depth of the crust and sensitive clay (Fig. 20b).F 5 has a shallower depth, thus it does not cross the slope base (Fig. 20c), but it still only originates from the middle part.In a progressive flow slide failure mechanism characterised by higher sensitivity (S t ) and a steeper strength gradient (k), only a relatively small portion of the entire failure process duration is attributed to F 1 .
In comparison to cases with lower strength gradients, a higher strength gradient leads to the formation of more conspicuous inclined shear bands within the crust layer.When the failure surface extends down the slope and beneath the crust layer, indicating a spread failure, the soils above the failure surface undergo movement and dislocation, and are squeezed due to interaction with the front and rear soils.This interaction results in the development of inclined shear bands between soil blocks situated above the failure surface, as described by Locat et al. (2011).In this study, two types of distributions for inclined shear bands are identified: (1) Type I, characterised by vague densely-distributed shear bands; and (2) Type II, characterised by loosely-distributed but clear shear bands, which form triangular shapes.However, due to the insensitive nature of the crust layer, there is no strain-softening, which means that even strong inclined shear bands do not exhibit apparent dislocation.Nevertheless, triangular soil blocks can form and create ploughing shapes, as noted in the work of Puzrin (2016).Furthermore, in Type II (e.g.Fig. 20d), some folds appear on the ground surface due to the presence of strongly inclined shear bands, a phenomenon observed in the northern California continental margin (Gardner et al., 1999).In contrast, Type I distribution presents a smoother ground surface (e.g.Fig. 19c).
Fig. 21 illustrates the detailed forming process of shear bands in Type II at t ¼ 10 s in Case 16 (Fig. 20b).The forming process for each triangular shape is consistent, and the sizes are similar.Unlike typical spreads, horsts and grabens do not arise from substantial dislocation in the crust layer.Nevertheless, the angle of the failure surface forming each triangle with its tip upwards (e.g. the acute   angle between shear bands 2 and 3 in Fig. 21a) approximates 45 , similar to the reported angle of a horst (Locat et al., 2011).After shear band 1 appears from the second flow slide (F 2 ), the bottom shear band 2 consistently emerges first, followed by shear band 3 (Fig. 21a).This process repeats to create the next triangle, and shear band 1 in the next triangle appears almost simultaneously with shear band 1 in the previous triangle.Through the velocity (v) mapping shown in Fig. 21b, inclined shear bands cut the crust into triangular blocks.Each triangular block has its own speed and is initially separated from other blocks.However, after some time, the velocities of these blocks converge, and they move as an integrated unit.It is worth noting that shear band 3 may not always become visible when shear band 1 appears quickly with closer spacing.In such cases, inclined shear bands exhibit Type I distribution.Type I distribution typically follows the end of Type II distribution, which occurs as the spread failure weakens and the collapse velocity decreases (Fig. 21a).Another scenario is when the collapse process of a spread failure is relatively slow and lasting, such as Case 14 in the case of low sensitivity (Fig. 19b and c), and Case 17 in the case of a low strain-softening rate (Fig. 22a), or when the failure is severe with the dominance of progressive flow slides, as in Cases 10, 11, and 4 (Fig. 15), resulting in only a Type I distribution due to spreads in the crust layer.

Effects of the rate of strain-softening
Cases 17, 15, and 18 in Table 1 depict three simulations with varying rate parameters of strain-softening (h) set at 5, 10, and 15, respectively.Fig. 22 displays the final deposit in plastic mapping format, while Fig. 23 shows the mapping of total displacement during the entire collapse for the three cases.As depicted in Fig. 22, a higher strain-softening rate leads to a longer diffused distance (L D ).This occurs because a larger h results in less unloading, which in turn produces a larger diffused distance, as reported by Locat et al. (2013).It is noteworthy that cases with faster strainsoftening, indicated by a larger h, are more prone to progressive flow slides, while cases with slower strain-softening are dominated by spread failure.In Case 17, only one spread failure, S 1 , originates from the initial flow slide, F 1 (Fig. 23aec), whereas progressive flow slides are evident in the other two cases.In Cases 15 and 18, F 2 appears earlier due to a higher strain-softening rate (Fig. 23e and g).F 3 in Cases 15 and 18 initiates from a closed position, but in Case 15, it does not extend across the slope bottom (Fig. 23f), and its horizontal span is smaller than F 3 in Case 18 (Fig. 23i).The simulated results indicate that a Type I distribution of inclined shear bands is more likely to appear in situations with a low strain-softening rate  as observed in Case 17 (Fig. 22a).Conversely, a Type II distribution is more probable in scenarios with a high strain-softening rate, as evident in Case 15 (Fig. 22b) and Case 18 (Fig. 22c).

Conclusions
The numerical simulations conducted using the N-PFEM method in this study have provided insights into post-failure collapses in sensitive clay on a long natural slope.The simulated results have revealed a variety of downward progressive landslide mechanisms, including single flow slides, translational progressive landslides, progressive flow slides, spread failure, or a combination of these mechanisms.The primary focus of this study has been to investigate the evolution of the downward progressive failure process and understand the effects of key parameters, namely sensitivity clay thickness (H s ), strength gradient (k), sensitivity (S t ), and strain-softening rate parameter (h).Several important conclusions have emerged from this investigation: (1) Four common failure mechanisms have been identified in this study, encompassing single flow slides, progressive flow slides, translational progressive landslides, and spread failures.These mechanisms often occur in combination during downward progressive landslides in sensitivity clays on extended natural slopes.(2) A local shear band, denoted as F 1 , consistently initiates near the embankment at the beginning of a collapse event.
Remarkably, its size and shape remain largely unaffected by the variations in H s , S t , and h.However, the strength gradient (k) does affects the shape of F 1 with higher k values resulting in a reduction in the size of F 1 .(3) Translational progressive landslides are particularly noteworthy when either all sensitive clays are located right below the crust layer or there is a very thin sensitive clay layer at the base.The former, characterised by high local and propagation speeds, poses significant hazards due to the instantaneous formation of a long failure surface.In the latter scenario, while the speed remains relatively low, the extensive range of diffused areas and the prolonged stabilisation time should not be underestimated as indicated in Case 2. (4) Flow slides are a common feature when all sensitive clays are situated right beneath the crust layer, often leading to spread failures underneath the crust.These failures are relatively perilous, given their rapid propagation and wide destruction area.
(5) In a failure which includes both spread and at least one flow slide, higher sensitivity, higher strain-softening rate, or lower strength gradient lead to a dominance of flow slides, whereas spreads are more pronounced.(6) Irrespective of whether all sensitive clays are beneath the crust layer or there is a very thin sensitive clay layer at the base, higher sensitivity, a greater strain-softening rate, and a lower strength gradient can all lead to increased final diffused distances.However, it should be noted that in the latter scenario, a high strength gradient may inhibit failure initiation due to the elevated strength of sensitive clays at the bottom.(7) Two types of inclined shear bands distributions are observed in this study: Type I, with multiple vague inclined shear bands; and Type II, with clear inclined shear bands forming triangular shapes.Type I usually appears in the crust layer by spreads due to low sensitivity or low strain-softening rate, it also could appear in a transient spread generated by a flow slide.Type I is also observed in a severe translational progressive landslide, penetrating the entire depth of the slope, such as Case 12. Type II usually appears in a faster and longer process spread failure in a case of higher sensitivity or strainsoftening rate, and it generally does not appear independently since its end will be distributed in Type I.
Some similarities between upward and downward progressive landslides have been identified, particularly regarding the effects of sensitivity, strain-softening rate, and strength gradient on final diffused areas.Notably, high sensitivity tends to correlate with an increased likelihood of flow slides, a trend supported by both numerical simulations and field observations.However, it is crucial to underscore the ongoing challenge in comprehensively studying the failure mechanisms of downward progressive landslides due to the limited research in this domain.While this study offers initial insights, it is essential to recognise that various other factors, such as soil properties, slope geometry, load distribution, and the positioning and extent of the sensitive clay layer, may also exert significant influence.Hence, it is imperative to give these factors further consideration in future investigations.

Fig. 4 .
Fig. 4. The basic steps of the N-PFEM: (a) Particle cluster; (b) Meshes of the computational domain; (c) Cells constructed through meshes; and (d) Particle cluster for the next step.
) "*" represents the dominant failure model in cases where multiple failure mechanisms occur.(3) Cases 12, 13, and 16 are not simulated until they reach stability due to boundary limits.(4) Type I e Multiple vague inclined shear bands; Type II e Clear inclined shear bands forming triangular shapes.(5)L D and U E represent the diffused distance and embankment displacement, respectively, as shown in Fig.6.It should be noted that both are considered in the horizontal direction (x-axis).

Fig. 5 .Fig. 6 .
Fig. 5. Evolution of undrained shear strength (su) in the numerical strain-softening model for different values of parameter h.

Fig. 12 .
Fig. 12. Final L D and U E against St in the case of k ¼ 1.2, h ¼ 10, and Hs ¼ 2 m.SF denotes the single flow slide, and TR denotes the translational progressive landslide."*" represents the dominant failure mode in cases where multiple failure mechanisms occur.

Fig. 13 .
Fig. 13.Evaluation of L D and U E against time in the case of k ¼ 1.2, St ¼ 2, and Hs ¼ 2 m.

Fig. 14 .
Fig. 14.Effects of strength gradient (k, in kPa/m) on the first shear band: (a) Sketch of the first shear band; and (b) Dimensions of the first shear band.

Table 1
Summary of typical cases in this study: Model setup and post-failure results.