Introduction

A retrogressive landslide is a type of geohazards that usually occurs in sensitive clays. A distinctive feature of this kind of landslides lies in the observed fact that an initial small local slope failure may lead to a catastrophic landslide proceeding rapidly in a retrogressive manner and eventually eroding seemingly stable areas that are kilometres away from the location of the initial failure. Representative retrogressive landslides include multiple retrogressive landslides (also called earthflows), translational progressive landslides and the spreads (Locat et al. 2011a; Cruden and Varnes 1996).

A multiple retrogressive landslide, as its name implies, involves a series of failures occurring sequentially (Locat et al. 2011a; Mitchell and Markell 1974; Quinn et al. 2011). Specifically, the landslide develops in a manner that the migration of mass from the initial slide results in a new backscarp which may be unstable as well. Failures of the newly generated backscarps occur until a final stable scarp is formed. The translational progressive landslide develops in a very different way; that is, a shear surface parallel to the ground surface develops first and then the soil mass above the formed failure surface moves downhill (Bernander et al. 2016; Wang 2019). In a spread, a shear surface along the ground surface forms as that in a translational progressive landslide. However, instead of moving downhill, the soil mass above the failure surface undergoes extension and dislocation producing blocks of intact clay called horsts and grabens (Locat et al. 2013; Geertsema et al. 2018; Locat et al. 2014).

Although computer simulation has been regarded as a powerful tool for slope stability analysis, the numerical modelling of a retrogressive landslide is still a challenging task. The conventional numerical approaches for slope stability analyses such as the limit equilibrium method and the finite element strength reduction method fail because of two issues. First of all, these approaches concern only the initiation of the landslide, for example the failure of the original slope, whereas for a retrogressive landslide, specifically for an earthflow, a series of collapses are triggered after the initiation of the landslide, which results in a large retrogression distance that neither the limit equilibrium method nor the finite element strength reduction method can predict. In addition, the strain softening feature of sensitive clay, which plays a key role in the forming of shear bands and consequently in the failure mode, cannot be captured by the mentioned conventional methods for slope stability analysis because of their assumption of perfect plastic soil behaviour. Indeed, a reliable numerical prediction of a retrogressive landslide and its consequences requires not only a sophisticated constitutive model for representing the complex behaviour of sensitive clay but also a robust numerical approach capable of simulating the entire process where large soil deformations occur and new free surfaces (e.g. new backscarps) emerge.

So far, numerous efforts have been devoted to reproducing the complete process of a retrogressive landslide. Recent contributions in this regard were made in (Dey et al. 2016; Wang et al. 2016; Dey et al. 2015) where the coupled Eulerian-Lagrangian (CEL) method available at the commercial software Abaqus and the material point method (MPM) were adopted for tackling the induced large deformations and the emergence of new free surfaces in the post-failure stage. In these cases, a rate-independent model in the form of the von Mises yield criterion with strain softening was applied to capture the progressive failure behaviour of sensitive clay. Instead, the particle finite element method (PFEM) that is a combination of the particle approach and the finite element method was employed to analyse retrogressive landslides in (Zhang et al. 2017) by adopting a rate-dependent elastoviscoplastic model constructed as a mixture of the Tresca model (with strain softening) and the Bingham model. More specifically, the part of Tresca yield criterion with strain softening in the model was used to capture the solid-like behaviour of sensitive clays. On the other hand, the Bingham plastic rheology was used since it is able to mimic the fluid-like behaviour (e.g. non-Newtonian flow) which is a typical feature of high-sensitivity clays, such as quick clays, when they are fully remoulded. Some key characteristics in retrogressive landslides have been reproduced successfully in those works (Dey et al. 2016; Wang et al. 2016; Dey et al. 2015; Zhang et al. 2017) including the development of progressive failure, multiple collapses, formation of horsts and grabens, etc. Despite that, these works focused on the qualitative studies of the mechanism of retrogressive landslides based on virtual landslide models. To the authors’ best knowledge, there have been so far very few quantitative studies of retrogressive landslides focusing on the comparison of numerical results with physical model tests or with real landslide cases in sensitive clay. In this paper, a case study on a retrogressive landslide that occurred in 2010 in Saint-Jude, Quebec, Canada, is carried out numerically based on the computational framework of the PFEM associated with an elastoviscoplastic constitutive model with strain softening.

The idea of the PFEM is that the Lagrangian finite element method is adopted for incremental analyses at each time step whereas the element nodes, after each incremental step, are considered free particles that may move freely and even separate from the domain they originally belong to. Computational domains and the associated meshes are constructed based on the particles when material undergoes large deformations. In other words, when modelling landslide problems, the PFEM is able not only to capture the free surface evolution (including the sliding mass movement and the emergence of new backscarps after slope failure) forthrightly owing to the Lagrangian description but also to circumvent mesh distortion issues that regard particularly strong shear deformation zones. The objective of this work is twofold: (a) to assess the performance of the computational framework in reproducing a real retrogressive landslide through detailed comparison between simulation results and data from site investigations and (b) to further analyse the mechanism and evolution of the Saint-Jude landslide through a numerical simulation with the eventual goal to improve the understanding of real retrogressive landslides.

The Saint-Jude landslide

The Saint-Jude landslide took place in a deposit of sensitive Champlain Sea clay (Lafleur and Lefebvre 1980) along the Salvail River on 10 May 2010 (Fig. 1). Leading to death of four people, the landslide had a retrogression distance of 80 m (defined as the distance from the initial crest of the slope to the final stable backscarp) and affected a section about 275 m long, parallel to the Salvail River and a total of 520,000 m3 soil mass (Locat et al. 2017). Both on-site and lab testing was carried out to obtain information on the stratigraphy and the geotechnical properties of the clays. The topographies of the slope before and after failure are taken from the digital elevation model built from aerial photographs. The chosen cross-section C-C′ before and after failure is illustrated in Fig. 2 for comparison.

Fig. 1
figure 1

The Saint-Jude landslide in 2010 (after (Locat et al. 2017)): (a) its location with dark grey area illustrating the Champlain Sea deposit and (b) aerial view. Red arrows represent the sliding direction

Full size image
Fig. 2
figure 2

Cross-section C-C′ of the slope profile: (a) three times vertical exaggeration, (b) before the Saint-Jude landslide and (c) after the Saint-Jude landslide. The triangles represent the identified horsts with C-i being the ith horst for cross-section C-C′ (modified from (Locat et al. 2017))

Full size image

The exact cause of the Saint-Jude landslide is unknown. Slope stability analyses were performed in (Locat et al. 2017; Locat et al. 2011b) along the transect C-C′. In (Locat et al. 2011b), it was shown that the initial failure based on drained analysis is very local whereas a large amount of clay would be disturbed if the slope failed under undrained conditions. These different failure modes lead to a substantial discrepancy in the factor of safety (FoS) determined under undrained and drained conditions. According to (Locat et al. 2017; Locat et al. 2011b), the analysis under undrained conditions using the limit equilibrium method provides the FoS of 2.16, whereas the FoS from drained analysis is close to one (more precisely 0.99 with the Bishop method and 1.03 with the Morgenstern-Price method). For detailed slope stability analyses under drained and undrained conditions including the used material parameters, we refer the reader to (Locat et al. 2011b). It was believed that the slope could have been initiated by small-magnitude triggering events such as the river erosion and the high artesian pore pressure under the river. For instance, one scenario could be that the river erosion may remove the materials at the slope toe leading to an immediate increase of the total stress at the potential shear band zone in the slope, while the high artesian pore pressure may increase the pore water pressure in the slope. Both of these two mechanisms may result in a decrease of FoS. The critical slip surface for the initial slope determined based on the drained analysis is illustrated in Fig. 2(a) (see the dashed blue curve).

The failure surface of the landslide was determined using cone penetration tests with pore water pressure measurement (CPTU test) in a series of locations (see for example the black dashed lines in Fig. 2(b), (c)). More specifically, from the CPTU profiles, one can derive the depth of the failure surface at each location (see e.g. black points in Fig. 2(a)) and then build the entire failure surface by a suitable interpolation. The reader is referred to (Locat et al. 2017; Locat et al. 2011b) for more details on the obtained CPTU profiles and the procedure to determine failure surfaces. According to the report (Locat et al. 2017), the debris of the retrogressive landslide can be divided into four zones (Fig. 2). The first zone occupies about 22% of the total landslide area. Here the soil mass comes from the initial bed and banks of the Salvail River that were pushed onto the opposite side of the river. Settled on top of the initial position of the Salvail River, zone 2 is about 20% of the landslide area and less fissured than zone 1. The ground surface of this zone is roughly horizontal with trees slightly inclined towards the backscarp. Blocks such as grabens and horsts were observed in both zone 3 (24% of the landslide area) and zone 4 (34% of the landslide area). While stratigraphy in zone 3 indicates that blocks did not rotate during the landslide motion, horsts in zone 4 were inclined between 0 and 17° to the horizontal. The sides of clay ridges formed by horsts in zone 4 are inclined around 45 to 80°.

As indicated in (Locat et al. 2017), the limit equilibrium method cannot explain the entire event. A better explanation of this geohazard and the associated failure mechanism requires more advanced calculation approaches. In this paper, the PFEM and the elastoviscoplastic constitutive model with strain softening introduced in (Zhang et al. 2017) are adopted for the back analysis of the Saint-Jude landslide.

Problem description

The aerial view and on-site investigation showed that the final profiles of three sections A-A′, B-B′ and C-C′, all along the landslide direction, after the landslide are very similar to each other. Therefore, in our study, the landslide is considered a plane strain problem, which is in line with the stability analysis performed in (Locat et al. 2017; Locat et al. 2011b). The numerical model in this study is constructed according to the cross-section C-C′ (Figs. 1(b) and 2(b)) whose critical slip surface (i.e. the blue dashed curve in Fig. 2(a)) determined using the limit equilibrium method under drained conditions was given in (Locat et al. 2017).

Governing equations

As the landslide process lasted only a few minutes and the permeability of clay is very low (Locat et al. 2017), the simulation is conducted by assuming undrained conditions, which is not in contrast with the fact that drained conditions were invoked to explain the slope instability, as discussed later on. The following governing equations have to be fulfilled on the volume Ω and the boundaries Γ:

  1. (a)

    The momentum conservation equation

$$ {\nabla}^{\mathrm{T}}\boldsymbol{\upsigma} +\mathbf{b}=\rho \overset{..}{\mathbf{u}}\kern1em \mathrm{in}\kern0.33em \varOmega $$
(1)
  1. (b)

    The elastoviscoplastic constitutive model with softening

$$ \left\{\begin{array}{l}\boldsymbol{\upsigma} =\boldsymbol{\uptau} +\eta {\dot{\boldsymbol{\upvarepsilon}}}^{\mathrm{vp}}\\ {}{\dot{\boldsymbol{\upvarepsilon}}}^{\mathrm{e}}=\mathbb{C}\dot{\boldsymbol{\uptau}};\kern1em {\dot{\boldsymbol{\upvarepsilon}}}^{\mathrm{vp}}=\dot{\lambda}\frac{\partial F\left(\boldsymbol{\uptau}, \kappa \right)}{\partial \boldsymbol{\uptau}}\\ {}\dot{\lambda}F\left(\boldsymbol{\uptau}, \kappa \right)=0;\kern0.66em \dot{\lambda}\ge 0;\kern0.66em F\left(\boldsymbol{\uptau}, \kappa \right)\le 0\end{array}\right. $$
(2)
  1. (c)

    The boundary conditions

$$ \left\{\begin{array}{l}\mathbf{u}=\overline{\mathbf{u}}\kern2.12em \mathrm{on}\kern0.66em {\varGamma}_{\mathrm{u}}\\ {}{\mathbf{N}}^{\mathrm{T}}\boldsymbol{\upsigma} =\overline{\mathbf{t}}\kern1em \mathrm{on}\kern0.66em {\varGamma}_{\mathrm{t}}\end{array}\right. $$
(3)

where

σ:

is the Cauchy stress

τ = σ:

if F(σ, κ) ≤ 0

b:

is the body force of clay

ρ:

is the density of clay

u:

is the displacement field

εe and εvp:

are the elastic strain and the viscoplastic strain, respectively, providing the total strain ε = εe + εvp

ℂ:

is the elastic compliance matrix

η:

is the viscosity coefficient

λ:

is the non-negative plastic multiplier

κ:

is the equivalent deviatoric plastic strain

\( \overline{\mathbf{u}} \) and \( \overline{\mathbf{t}} \):

are the prescribed displacements and external tractions on boundaries

N:

consists of components of the outward normal to the boundary Γt

∇T:

is the transposed gradient operator.

The superposed dot represents differentiation with respect to time.

Following (Zhang et al. 2017), the Tresca yield criterion

$$ F\left(\boldsymbol{\upsigma}, \kappa \right)=\sqrt{{\left({\sigma}_{xx}-{\sigma}_{yy}\right)}^2+4{\sigma}_{xy}^2}-2{c}_{\mathrm{u}}\left(\kappa \right) $$
(4)

is used in which the strain softening behaviour of sensitive clay is taken into account by reducing the undrained shear strength with the equivalent deviatoric plastic strain based on a two-segment linear function as shown in Fig. 3. The equivalent deviatoric plastic strain is defined as \( \kappa =\sqrt{0.5{e}_{ij}^{\mathrm{p}}{e}_{ij}^{\mathrm{p}}} \) according to (Troncone 2005) with \( {e}_{ij}^{\mathrm{p}} \) being the deviatoric plastic strain tensor. The present elastoviscoplastic model reduces to the classical rate-independent elastoplastic model in the case of η = 0. The interested readers are referred to (Zhang et al. 2017) for more details on the abovementioned governing equations.

Fig. 3
figure 3

Variation of the undrained shear strength (cu) with equivalent deviatoric plastic strain (κ)

Full size image

Clay properties

According to (Locat et al. 2017; Locat et al. 2011b), soils at the elevation less than 2 m (see the vertical axis of Fig. 2(b)) are much stiffer than the soils above. Thereby, only soil layers of elevation more than or equal to 2 m are considered in our simulation, including the top crust (a dense sandy layer) and the underneath deposit of firm grey sensitive clay (very uniform with some silt). In the simulation, the crust is considered non-sensitive with shear strength of 75 kPa which is in line with the parameter used for undrained analysis in (Locat et al. 2011b). Field vane testing (Locat et al. 2011b) showed that the intact shear strength increases linearly with the depth of the deposit from 25 to 65 kPa. Thus, the peak shear strength used for the sensitive clay in our simulation is the average value (45 kPa) with the residual shear strength being 1.6 kPa, which is the strength of the clay near the bottom of the deposit when it is fully remoulded. According to (Locat et al. 2017; Locat et al. 2011b), the unit weights of the crust and the sensitive clay are 18.6 and 16 kN/m3, respectively, and the Poisson’s ratio for both the crust and the sensitive clay is 0.49 to constrain the volume change (undrained conditions). The values of the reference equivalent deviatoric plastic strain and the viscosity coefficient were not reported in (Locat et al. 2017; Locat et al. 2011b). In our simulation, the reference equivalent deviatoric plastic strain is \( \overline{\kappa}=1.6 \) calculated from the lab testing data for Canadian sensitive clay reported in (Quinn et al. 2011). The viscosity coefficient for the crust and the deposit of sensitive clay is 1000 Pa ⋅ s which is within the suggested range 100–1499 Pa ⋅ s for case studies of landslides (Edgers and Karlsrud 1982; Johnson and Rodine 1984) and in line with the viscosity for reproducing the Storegga slide in a marine sensitive clay (Gauer et al. 2005).

Solution scheme—PFEM

A variety of numerical techniques are available for solving the above governing equations. In this work, the set of equations is solved using the finite element method (FEM) in mathematical programming as suggested in (Zhang et al. 2017). In this computational framework, the governing equations after finite element discretization are reformulated as an equivalent second-order cone programming problem, based on the Hellinger-Reissner variational principle, which is then solved by employing an advanced optimization algorithm—the primal-dual interior point method. For more details of the solution scheme, readers are referred to (Zhang et al. 2017) in which the correctness and robustness of the solution scheme have been illustrated.

Since the used FEM is based on Lagrangian description, it fails in modelling the entire failure process of the retrogressive landslide. This is so because the clay undergoes very large deformations that result in severe mesh distortion if Lagrangian FEM is used. Moreover, the change of geometry may lead to a complex free surface evolution, for example to the generation of new backscarps in the collapse, which cannot be captured by the traditional Lagrangian FEM. To tackle these issues, the version of the PFEM introduced in (Zhang et al. 2017) is also adopted in this work. The PFEM is a mixture of the Lagrangian FEM and the particle approach. The basic idea behind the PFEM is that mesh nodes are regarded as free particles that can move freely while governing equations are still discretized and solved using Lagrangian FEM incrementally (Zhang et al. 2017; Oñate et al. 2004). For the concerned problem, the basic analysis steps for the PFEM in a typical time interval [tn,  tn + 1] are as follows (see also Fig. 4):

  1. (a)

    Solve the governing equations on mesh Mn, using the Lagrangian FEM;

  2. (b)

    Update the mesh nodes to their new position according to the computed incremental displacement and erase the mesh topology leaving behind a cloud of free particles Cn + 1;

  3. (c)

    Identify the computational domain Ωn + 1 based on Cn + 1 using the alpha-shape method;

  4. (d)

    Discretise the domain Ωn + 1 to obtain the new mesh Mn + 1;

  5. (e)

    Map the state variables, such as velocities, accelerations, stresses and strains, from the old mesh to the new mesh and then step to the next time interval.

Fig. 4
figure 4

Illustration of basic steps of the PFEM in a typical time increment

Full size image

A key feature of the PFEM lies in its inheritance not only of the capability of a particle approach for handling general large deformation problems but also of the solid mathematical foundation of the FEM. To date, the PFEM has been used successfully for modelling landslides and their relevant large deformation problems such as the granular flow (Dávalos et al. 2015; Zhang et al. 2014; Cante et al. 2014; Zhang et al. 2016), the debris flow and its interaction with barriers (Franci and Zhang 2018), a flow-like landslide (Zhang et al. 2015), landslide-generated waves (Salazar et al. 2016; Cremonesi et al. 2011), fluid-soil-structure interactions (Oñate et al. 2011), submarine landslides and their impact to nearby subsea pipelines (Zhang et al. 2019), just to mention a few. The possibility of the used PFEM for modelling retrogressive landslides has been explored in (Zhang et al. 2017). In this work, it is used to study the 2010 Saint-Jude retrogressive landslide qualitatively and quantitatively.

Simulation results and discussions

To trigger the landslide, the clay enveloped by the critical slip surface determined from the drained slope stability analysis performed in (Locat et al. 2017; Locat et al. 2011b) is set to be fully remoulded. This critical slip surface is selected because the FoS from drained analysis (implying that the excess pore water pressure is zero) is much lower than that from the undrained analysis (where the clay is impermeable and the excess pore water pressure cannot be dissipated). However, after the landslide is triggered, the time interval from the mass transport to the final deposition is very short, and therefore, given that the clay permeability is very low, the effect of water seepage in clay can be neglected. This explains why the post-failure analysis in this work is modelled under undrained conditions. In our simulation, the profile of the deposit is set to be sufficiently long to remove the boundary effect. The bottom of the sensitive clay is fully fixed. A total of 17,773 mixed triangular finite elements (36,224 mesh nodes) are used to discretize the initial domain of the deposit. Mixed elements have been introduced in (Zhang et al. 2017), and the term “mixed” means that in each element displacements and stresses are treated as independent fields: Quadratic shape functions are used for displacements whereas linear shape functions are used for stresses. The time step adopted in the simulation is Δt = 0.05  s.

Failure evolution and its implication

The complete failure process of the Saint-Jude landslide from our simulation is shown in Fig. 5 where the contour represents the distribution of the equivalent deviatoric plastic strain κ. Overall, the Saint-Jude landslide is in a form of a multiple rotational retrogressive landslide according to the simulation results. The mass involved in the failure of the initial slope flows into the river leading to the yielding of the clay at the bottom of the newly generated backscarp (Fig. 5(a)). The accumulation of plastic deformation in the backscarp at this stage is relevantly slow since the clay from the first failure is restored in the river that is just in front of the new backscarp. Nevertheless, the failure band still develops gradually from the bottom of the deposit to its top surface resulting in the first retrogressive failure (RF) as shown in Fig. 5(b). A relatively large amount of clay, more than twice as much as that involved in the failure of the initial slope, is involved in this failure. The gravitational potential energy of the clay is transformed into its kinematic energy, and thus, the clay from the first RF pushes the remoulded clay in front to further move forward. An apparent strength reduction is observed for the clay in the failure band while the clay outside the failure band is almost intact and migrates with little deformation. The migration of the clay from the first failure results in a crater in front of the backscarp (Fig. 5(c)), and similarly, the second RF is then triggered and behaves as a rotational failure at the early stage (Fig. 5(c)). The amount of clay disturbed by the second RF is even larger. As can be seen in Fig. 5(c), the shear band of the failure at this time is much longer, which is, to a large extent, due to the increase of the deposit thickness at the distance of 140 m. During the migration of the clay involved in the second RF, the stress intensity at the position where thickness increases evokes two additional shear bands propagating downwards making the second RF evolve from a rotational collapse to a spread process (Fig. 5(d), (e)). The two additional shear bands divide the block in zone 3 into grabens and horsts (see the triangles in Fig. 5(f)). These simulated phenomena have two implications: (1) the evolution of a rotational failure into a spread is possible in a thick sensitive clay layer during the migration process and (2) a retrogressive landslide can be the combination of a multiple retrogressive failure and a spread. While the clay from the second RF migrates, the third RF occurs in the mode of a rotational failure (Fig. 5(e)) resulting in a final stable backscarp shown in Fig. 5(f). Notably, the block involved in the third failure also breaks due to the new emerging shear band, resulting in two parts in the shape of a graben and a horst. It is thus concluded that not only a spread but also a rotational retrogressive failure can give rise to grabens and horsts. Figure 5 (f) shows that grabens and horsts can be observed in both zone 3 and zone 4, which agrees with the observations reported in (Locat et al. 2017).

Fig. 5
figure 5

Contours of the equivalent deviatoric plastic strain, at (a) t = 4 s, (b) t = 56 s, (c) t = 65 s, (d) t = 69 s, (e) t = 76 s and (f) t = 87 s. RF denotes a retrogressive failure. The grey solid line indicates the initial profile of the slope, while the red dashed line and the black circles in (f) correspond to the final deposit profile and the failure surface determined by on-site tests (Locat et al. 2017), respectively. The undrained strength is at the peak value if κ = 0 (the materials in blue) while it is at the residual value when κ ≥ 1.6 (the materials in red)

Full size image

Site investigations also show that blocks in zone 3 did not rotate during the landslide whereas horsts in zone 4 were inclined between 0 and 17° with respect to the horizontal. These observations can be explained by our simulation where the migration of the clay in zone 3 is in the form of a spread where horsts usually undergo very limited rotation, whereas zone 4 results from a typical rotational retrogressive failure where blocks rotate considerably during the migration.

The final deposit profile of the Saint-Jude landslide (red dashed line) and the failure surface of the landslide determined by on-site testing (black circles) are also depicted in Fig. 5(f) for comparison. As can be seen, the final run-out distance and the retrogression distance from our simulation are 55 m and 145 m (80 m if measured from the crest of the initial slope), respectively, which agrees well with the observation data reported in (Locat et al. 2017). A satisfactory agreement in terms of the shape of the final deposit is also achieved. Simulation results show that the entire landslide area is formed by four zones, and this is consistent with the observations reported in (Locat et al. 2017). Moreover, our simulation results indicate that a single zone of the landslide area is originated from an individual failure. More precisely, the first zone is from the initial failure, and the second zone is from the first RF and so on. Further, the failure surface from our simulation matches that one determined from the piezocone penetration tests (CPTu). In conclusion, the qualitative and quantitative agreements validate the adopted computational framework for modelling real retrogressive landslides in sensitive clay.

Failure mechanism

Shear stress fields τxy in the deposit at different instants are illustrated in Fig. 6. These time instants refer to the moment when a new retrogressive failure (see the red curves) is triggered. The black curves in Fig. 6 represent the shear stress τxy at the bottom of the deposit (namely at the elevation of 2 m). It is worth noting that, since the Tresca yield function is used, both the norm of shear stress τxy and the norm of (σxx − σyy) contribute to the yielding of the clay. Nevertheless, only shear stress field is discussed here since our simulation results show that the magnitude of |σxx − σyy| is much lower than the magnitude of |τxy|.

Fig. 6
figure 6

Contours of shear stress at (a) t = 4 s, (b) t = 56 s, (c) t = 65 s, (d) t = 76 s and (e) t = 87 s. Black curves correspond to the shear stress on the bottom of the deposit, and red curves represent the failure surface of the new retrogressive failure (RF)

Full size image

As seen in Fig. 6(a), when a new backscarp (the red curve in Fig. 6(a)) is created owing to the removal of clays involved in the failure of the initial slope, the shear stress at the bottom of the sensitive clays, represented by the solid black curve, behind the new scarp increases to a very high value. The clay there yields with its strength being further reduced, which results in the emergence of a new shear band. The shear band propagates from the bottom of the clay layer to the top surface of the deposit, which triggers the first retrogressive failure (RF) as shown in Fig. 6(b). The second and third sequential collapses (Fig. 6(c), (d)) follow the same failure mechanism.

When a new RF is triggered, the shear stress on the bottom of the deposit on the left of the new RF surface is very low because the clay there has undergone large plastic deformation and consequently possesses a very low shear strength. On the other hand, the shear stress on the right-hand side of the new RF surface is much higher (slightly below 45 kPa). This is because the clay in this area is slightly disturbed at the moment. The shear strength of the clay is still high, and thus, it can sustain a relevantly large shear stress.

It is also notable that the curves of the shear stress on the bottom of the clay deposit fluctuate somewhat in the area in front of the new RF surface (e.g. at the distance between 50 and 150 m at t = 65 s as shown in Fig. 6(c), particularly at the distance of 70 m and 100 m). This fluctuation stems partly from dynamic effects. But partly it is due to the existence of some intact clay or partially remoulded clay in these areas which can still sustain a high shear stress. Indeed, partially remoulded clay is observed at the distance around 70 m and 100 m at t = 65 s in Fig. 5(c). As the result of the continuous migration of the clay, the strength of the clay further reduces, and thus, the fluctuation of the shear stress in these areas is much less, for example at t = 76 s and 87 s as shown in Fig. 6 (d) and (e), respectively.

Kinematics

The velocity of clay in the landslide process is depicted in Fig. 7 where the meshes are also plotted. As seen, meshes of good quality are maintained in the simulation. The moving direction of the clay is generally parallel to the basal failure surface of the landslide, given by the solid black curves.

Fig. 7
figure 7

Velocities of clays in m/s at (a) t = 4 s, (b) t = 56 s, (c) t = 65 s, (d) t = 69 s, (e) t = 70.5 s, (f) t = 76 s and (g) t = 87 s. Also the final configuration is depicted in (g)

Full size image

Due to the rotational RF, the gravitational potential energy of the evoked clay is transformed into its kinematic energy, and thus, a relatively high speed is observed for the clay involved in the new RF. For the Saint-Jude landslide, the mass involved in a former failure is almost static or possesses a very low migration speed when a new RF is triggered. Thus, the mass evoked by the new failure usually possesses the maximum migration speed and pushes the mass in front to move forward (Fig. 7(b), (c), (f)). However, a huge amount of mass is disturbed in the second RF, which releases a large amount of gravitational potential energy. The wave soon propagates to the leading front of the landslide causing the front to possess the maximum migration speed.

Additionally, as the mass from the first RF is still almost at t = 65 s (Fig. 7(c)), the block disturbed in the second RF somewhat ploughed into the front mass which can be also identified in Fig. 5(d) (at the distance of around 90 m). Moreover, some small undisturbed areas are observed below the new failure surface when the second and third RFs occur, for example area A in Fig. 7(c) and area B in Fig. 7(f), because these RFs at the beginning are of a rotational failure type. The clay velocity in these areas is zero when the failure is triggered. Soon after, area A is eroded owing to the migration of mass, whereas area B is maintained until the end of the landslide. This indirectly supports the previous explanation of the fluctuation of the bottom shear stress in front of the new generated backscarp.

Effects of viscosity for Saint-Jude landslide modelling

Since the clay viscosity was not calibrated in (Locat et al. 2017), the viscosity coefficient of the clay used in our simulations (1000 Pa ⋅ s) is determined by back analysis. As shown in Fig. 8(a), the landslide leads to a run-out distance of 55 m and a retrogression distance of 145 m (80 m if calculated from the crest of the initial slope) that agree well with the reported data if the adopted viscosity is 1000 Pa ⋅ s. To further investigate the influence of the viscosity on the Saint-Jude landslide, a parametric study was carried out with viscosity being 500 Pa ⋅ s, 200 Pa ⋅ s and 0 Pa ⋅ s, respectively. The final deposits of the Saint-Jude landslide from the simulations are illustrated in Fig. 8. It is shown that decreasing the viscosity results in an increase of both the run-out distance and the retrogression distance. When the viscosity is neglected (Fig. 8(d)), the final run-out distance from the simulation is 75 m (36.4% greater than the reported one) and the retrogression distance is 170 m (17.2% greater than the reported one). When the viscosity equals 500 Pa ⋅ s or 200 Pa ⋅ s (Fig. 8(b), (c)), the final depositions obtained from the simulations are almost the same. The run-out distance is slightly larger than the reported data, and the retrogression distance matches the reported data although a shear band emerges at the location from 175 to 220 m which may lead to a further retrogression distance of 25 m (Fig. 8(b), (c)). In general, for the Saint-Jude landslide, the failure mode and the retrogression distance induced by each retrogressive failure are independent of the viscosity. It is thus concluded that the value of the clay viscosity ranging between 200 and 1000 Pa ⋅ s may reproduce a reasonable final deposition for the Saint-Jude landslide, which is a value within the range of viscosity for real landslide simulations (100–1499 Pa ⋅ s) suggested by Edgers and Karlsrud (Edgers and Karlsrud 1982) and Johnson and Rodine (Johnson and Rodine 1984). Nevertheless, it is noticeable that the used viscosity is much greater than the viscosity calibrated in lab testing for Canadian sensitive clay in (Locat and Demers 1988) that varies from 0.007 to 0.2 Pa ⋅ s. The discrepancy between the viscosity value determined in the lab testing and that used to predict a reasonable run-out distance for sensitive clay landslides is still an open question where there is no consensus since there is a variety of influencing factors such as the accuracy of the calibration method for viscosity and the used rheological models. To explain this discrepancy, research works on the synergistic investigations from both physical and numerical modelling are required.

Fig. 8
figure 8

Final depositions of the Saint-Jude landslide obtained from the PFEM modelling with clay viscosity being (a) 1000 Pa ⋅ s, (b) 500 Pa ⋅ s, (c) 200 Pa ⋅ s and (d) 0 Pa ⋅ s. Maps represent the equivalent deviatoric plastic strain

Full size image

Conclusions

In this paper, a case study on the Saint-Jude landslide that occurred in 2010 was carried out numerically using a new computational framework. The PFEM was adopted to tackle the issues resulting from large material deformation while the behaviour of sensitive clays was described by an elastoviscoplastic constitutive model that is a mixture of the Tresca model with strain softening and the Bingham model. The investigation of progressive failure in Saint-Jude landslide was carried out with emphasis on the failure mode and mechanism, the kinematics and the final deposition. Detailed comparisons between the modelling results and the available data were conducted. The main conclusions are drawn as follows:

  1. (1)

    The PFEM using the elastoviscoplastic model with strain softening is a valid computational framework for simulating real retrogressive landslides in sensitive clay. It has not only qualitatively reproduced the failure mode and the progressive failure process of the Saint-Jude landslide but also quantitatively estimated its final run-out distance, retrogression distance and failure surface. Thus, this computational framework can be used for hazard estimation of real landslides in sensitive clay.

  2. (2)

    The Saint-Jude landslide is a type of multiple retrogressive landslides (also called earthflows) where a total of three retrogressive failures (RFs) were triggered after the initial failure of the slope according to our simulations. The initial failure and the subsequent three RFs led the landslide into four zones which coincide with the observations.

  3. (3)

    Our simulation results indicated that a retrogressive landslide can be the combination of multiple types of retrogressive failure. This is supported by the evidence that the second rotational retrogressive failure in the Saint-Jude landslide evolves into a spread in the migration process.

  4. (4)

    Owing to the change of geometry, a spread may occur in a deposit of a thick layer of sensitive clays. This is demonstrated by the second retrogressive failure in the Saint-Jude landslide from our simulation. The stress intensity at the location where there is a sudden increase of clay thickness leads to additional shear bands which further decompose the moving block, make the failure evolve as a spread during migration and sequentially produce multiple grabens and horsts in the migrating clay.

  5. (5)

    Fluctuations of the shear stress on the bottom of sensitive clay will occur in front of the new failure surface because of the existence of intact or partially remoulded clay in these areas which may still sustain large shear stress. However, the fluctuation will be alleviated soon afterwards because the un-remoulded clay will be further eroded in the migration process.

  6. (6)

    The clay block stemming from a rotational failure may be further decomposed into two smaller parts exhibiting as a graben and a horst in the process of migration. In other words, a rotational retrogressive failure (not just the spread) in sensitive clay can produce grabens and horsts.

  7. (7)

    The horsts and grabens are blocks of intact clays that are usually under translation with very limited rotation in the process of migration if they are from a spread but rotate considerably if they derive from a rotational retrogressive failure.

  8. (8)

    In a multiple retrogressive landslide, usually the materials involved in a new retrogressive failure rather than those at the leading front possess the largest migration speed. However, if a huge amount of gravitational potential energy is released due to the new collapse, the wave may propagate fast to the leading front of the landslide and make it possess a very high migration speed.

  9. (9)

    Neglecting the viscosity of sensitive clays leads to an overestimation of both final run-out distance and retrogression distance for real retrogressive landslides. The modelling results for the Saint-Jude landslide suggest a value between 200 and 1000 Pa ⋅ s for producing reasonable run-out distance and retrogression distance. Although the suggested value is within the range advised by other researchers for real landslide modelling, it is much greater than the viscosity of Canadian sensitive clay determined by lab testing. Thus, synergistic investigations from physical and numerical modelling are required to provide a clear answer to this discrepancy issue.