A GIS framework for prediction of basin‐wide post‐failure geometry of confined subaqueous landslides

Stability analysis of submarine and sub‐lacustrine slopes requires special attention to the strain‐softening nature of clay‐rich sediments. Utilizing the often‐used infinite slope assumption using the peak strength is inadequate for the task at hand as the well‐known limit equilibrium method may provide non‐conservative estimate of marine slopes’ stability. In this paper, a GIS based framework is developed to account for a more suitable shear band propagation slope‐failure criterion for basin‐wide stability calculations. After determining the slab failure prone areas, the framework employs analytical criterion for retrogressive failure above the slab and depth‐integrated finite element method for calculating the size of confined slides at the toe. This allows predicting the landslide's post‐failure configuration and the overall landslide coverage zone. The framework is first benchmarked to prove its competence and then applied to the well documented Zinnen Slide and the southeastern banks of the Küssnacht Basin, Lake Lucerne, Switzerland. The results of the application show good agreement with the observed basin geomorphology.


INTRODUCTION
The evaluation of a slope susceptibility to sliding in submarine environment is a highly complicated task. The development of offshore infrastructure, given the high costs of deep-sea investigations, requires a solid prior knowledge of the risks involved in the preliminary design. Due to the large dimensions of submarine basins, site specific analysis is limited owing to the effort required to conduct such studies. Geographic Information Systems (GIS) are a well-established tool aimed at conducting studies over large areas, which has been used extensively to evaluate terrestrial, submarine, and sub-lacustrine slope stability, for example, Strasser et al., 1 Rushton et al., 2 Collico et al. 3 and many more. The ability of this tool to apply analytical criteria to slope failure using raster algebra has made it very attractive to use for analyzing large areas, both for its simplicity and time efficiency.
One of the most common methods to calculate the deterministic factor of safety (FoS) of slopes is the limit equilibrium method (LEM). Shallow-seated landslides' stability is often calculated assuming infinite slope conditions, 4,5 while stability of slides over deep-seated slip surfaces is estimated using various methods of slices. 6,7 In the LEM, linear momentum conservation between the driving forces (induced by gravity and/or seismic loading) and stabilizing forces (mobilized strength of the soil) is sought. If the ratio between the stabilizing forces and the driving forces (i.e., FoS) drops below unity, the slope is considered unstable. While in general quite powerful, the LEM reaches its limits in the case of strain-softening soils. Using the peak strength may provide non-conservative results due to overestimation of resistance. Considering the residual strength, on the other hand, may lead to over-conservative or non-plausible designs, as the residual strength of the soil may be lower than the shear stresses induced by gravity, leading the FoS to drop below unity even for existing in-situ conditions.
Weak layers interbedded within marine and lacustrine sediment horizons that exhibit strain-softening material behavior are known to serve as control for vast submarine slope failures. [8][9][10][11] Studies show how the presence of strain-softening weak layers overlain by layers of stronger sediments may provoke progressive growth of slip surfaces under certain conditions. [12][13][14][15][16][17][18][19] At the steepest parts of the slope, a certain length of the weak layer may soften if the driving shear stresses suddenly exceed the peak strength (for example due to seismic shaking or increase in pore water pressures), which Puzrin et al. 20 define as a primary failure. If the softened length exceeds a critical value, the slip surface grows catastrophically in an unstable fashion. 12,13,15 Ultimately, upslope active and downslope passive soil failure will take place and the slope will experience a global slab failure (secondary failure, after Puzrin et al. 20 ).
Puzrin et al. 15 and Zhang et al. 17 have derived failure criteria capable of accounting for nonlinear slope geometry, namely slopes which can be approximated by an exponential profile. Such approximation for submarine slope geometry was found to be in good agreement with a large number of seismic profiles of the world's continental margins. 21 Puzrin et al. 15 and Rushton et al. 2 have shown that application of the so called 'Shear band propagation (SBP) approach' results in lower FoS when compared to LEM using the peak strength, leading LEM to provide non-conservative higher safety margins. Moreover, the unstable area of the slope may be larger when considering the possible extent of the shear band propagating into a more stable ground.
Regardless of the method used to derive the spatial distribution of FoS, the area in which FoS equals or drops below unity marks only the initial dimensions of the slab failure. Upon continuous downward displacement of the sliding body, the stable top and bottom portions are subjected to increasing unloading and loading, respectively. If a certain drop of the ground surface is reached, tertiary failure mechanisms may be initiated characterized by recurring upslope retrogressive active failure blocks. Downslope, the landslide may either emerge and run-out or plough into the stable ground, creating recurring passive failure blocks until all energy is dissipated in plastic work. These tertiary failure mechanisms may increase the dimensions of the affected area beyond the initial unstable zones and their consideration is crucial for adopting appropriate design approach and mitigations. 22 Over the years, a considerable amount of work has been devoted to the understanding of these mechanisms and the estimation of their potential dimensions. Among the methods used were analytical energy-based criteria 20,23-26 and limit equilibrium [27][28][29] and various numerical methods such as large deformation finite element analysis, 28,30-32 material point method, 33,34 smooth particle hydrodynamics, 35 particle finite element method, 36 depth-integrated methods [37][38][39][40][41][42] and more.
Rushton et al. 2 introduced the criterion by Puzrin et al. 15 into GIS applications in order to conduct a geohazard study for the Azeri-Chirag-Gunashli development in the Caspian Sea. In their work however, several challenges have risen in the data analysis, which required making certain assumptions and manual intervention in the process. Most of these challenges originated from the analytical criterion being formulated for a 2D problem, that is, plane-strain conditions, and not real 3D bathymetries. In the presented here work, these challenges are resolved through implementation of few additional steps. Beyond failure initiation and shear band propagation, this paper presents a geomechanics based approach to estimate the landslide final dimensions. Firstly, the relevant theory required to estimate failure initiation dimensions of the failed slab is presented. Secondly, the analytical criterion of Zhang et al. 28 is implemented for determining the retrogression limit above this slab. Finally, a newly developed depth-integrated finite element method (DIFEM) with adaptive remeshing is developed to simulate the landslide evolution and to assess its downslope extent. After the methodology is presented, it is applied to an artificial slope of an exponential profile to test its functionality. After all the framework components are validated and benchmarked, the suggested framework is applied to a real-life bathymetry of Lake Lucerne, in central Switzerland, focusing on the well-documented case of the Zinnen Slide and the southeastern banks of the Küssnacht Basin.

The shear band propagation (SBP) approach
According to the limit equilibrium method, the FoS is defined by the ratio between the maximum shear resistance and the mobilized shear stresses where , and ℎ are the peak strength of the soil, the gravitational shear stress and a pseudo-static earthquake induced shear stress, respectively, given by where ℎ is the depth to the shear surface, ′ is the effective unit weight of the soil, is the slope inclination and ℎ is the pseudo-static earthquake coefficient is the peak ground acceleration, is the gravitational acceleration, is the total unit weight of the soil and is a pseudo-static seismic coefficient multiplier.
The peak undrained shear strength of marine sediments is often considered proportional to the consolidation pressure 43 by where is the undrained shear strength coefficient for normally consolidated sediments. The SBP approach uses the definition of the shear stress ratio, where and are the peak and residual strength of the weak layer, respectively, defined by with being the undrained shear strength increase coefficient of the weak layer and is the sensitivity of the weak layer.
The shear stress ratio allows distinguishing between three different zones of the slope, namely the unstable zone where the driving shear stresses exceed the peak strength (i.e., + ℎ ≥ or ≥ 1); the quasi-stable zone where the driving shear stresses are lower than the peak strength but exceed the residual ( ≤ + ℎ < or 0 ≤ < 1) and the stable zone where the driving shear stresses are lower than the residual ( + ℎ < or < 0). F I G U R E 1 (A) Exponential slope geometry and B) the different zones according to the SBP approach.
Puzrin et al. 15 investigated the process of SBP in submarine slopes approximated by an exponential profile of the form where and are the vertical and horizontal coordinates, respectively, is the maximal inclination of the slope and is half the total height drop ( Figure 1A). For derivation of the failure criterion, they considered an elastic sediment layer underlain by a strain-softening weak layer. A weak zone was assumed to form at the steepest part of the slope where the combined action of gravitational and earthquake-induced shear stresses overcome the peak strength, that is, in the unstable zone, where ≥ 1 ( Figure 1B). They showed that if the initial length of the unstable zone ( ) exceeds a critical length ( ), failure is initiated from the unstable zone and the shear band propagates catastrophically (under existing loads) to the borders between the stable and quasi-stable zones ( = 0), that is, over the entire length ( Figure 1A). The failure criterion is given by wherēand̄are the average shear stress ratio and inclination of the unstable zone, respectively. and are the planestrain loading and unloading moduli of the soil and̄is the characteristic displacement. The strength of marine sediments is often assumed to degrade linearly with the plastic deformation, , 44 in which casē= ∕2, where is the plastic shear deformation required to remold the soil from peak to residual strength. If the material exhibits clear nonlinear behavior, such as exponential strength softening, 45̄= 95 ∕3 where 95 is the value of plastic deformation required to cause 95% of the strength reduction. 46 In case ≥ the propagation of the shear band into the unstable zone will cause softening of the weak layer outside the failure initiation zone. Puzrin et al. 15 therefore suggested that the FoS over the entire length should be calculated using the residual strength, as opposed to LEM which utilizes the peak strength (Equation (1)) over the length of only. It follows that not only does the SBP approach yield a lower FoS comparing to LEM, but it also results in predicting a greater length of the initial sliding body.

2.2
The extent of retrogressive failure Zhang et al. 28 investigated the conditions for upslope retrogressive failure to take place by comparing the initial driving force (ˆ0), the initial weak layer resistance (ˆ0) and the initial sliding layer resistance (ˆ0). They distinguished between three scenarios, two of which are relevant to the present study, namely (a) high-strength sliding layer (ˆ0 ≥ 0 ) and (b) medium-strength sliding layer (ˆ0 ≥ˆ0). In the case of a sliding layer of low strength, the shear band cannot propagate into the stable zone and therefore considering the sliding layer to be of high strength provides a safer estimate of the extent of the failure zones. In both cases (a) and (b) the driving force is sufficiently large to cause initial failure, that isˆ0 ≥ˆ0 for case (a) andˆ0 ≥ˆ0 for case (b). The initial normalized sliding layer resistance is given byˆ0 where 0 is the coefficient of lateral earth pressure at rest, and are the characteristic length and stability number, respectively, calculated by wherēis the average peak strength of the sliding layer. If the maximal normalized driving force,ˆ, exceeds the initial normalized resistance of the sliding layer,ˆ0, a head scarp will form at the border between the quasi-stable and stable zones, that is = 0 (Figure 2A). Mathematically this is expressed as wherē= 0 is the average shear stress ratio in the interval ∈ {0, = 0 } and̂= 0 = =0 ∕ . The necessary condition for retrogressive failure 28 becomes with given by Equation (13). After initial failure, downward displacement of the failed slab causes formation of a graben at the position of the initial head scarp ( Figure 2B). Downward displacement of the graben and the subsidence of sediment thickness unloads the head scarp. Sufficient subsidence of the sediment thickness nucleates retrogressive failure by reducing the supporting stress to the active failure stresses. The critical sediment thickness ( Figure 2C) is given by where is the sensitivity of the sliding layer. Considering equilibrium of the failed blocks, the retrogression length, , can be calculated accordingly wherēis the current average shear stress ratio over the length ( Figure 2C). For simplicity, Zhang et al. 28 suggested that the current shear stress ratio distribution should be calculated using the critical height, ℎ Note that the calculation of the retrogression length involves an iterative process becausēis calculated over this length.

2.3
Downslope extent of failure using a mixed implicit-explicit dynamic depth-integrated finite element method (DIFEM)

2.3.1
Formulation of sliding layer bar elements The depth-integrated method is used to model terrestrial and submarine landslides in cases where the thickness of the sliding body is significantly smaller than its length. [37][38][39][40][41][42] The small thickness to length ratio of such slides allows the stresses and strains of the sliding body to be averaged over its thickness. As a result, the numerical model can be represented by a set of finite bar elements connected via mutual nodes ( Figure 3). This modeling technique is beneficial compared to  other numerical methods, such as the standard large deformation finite element or material point method, because of its short runtime and low computational resource requirements, which are advantageous for conducting studies that require a significant number of computations.
The bar elements are formulated such that each node has a single degree of freedom corresponding to the slope-parallel displacement, . The index i stands for the element/node number and the index j for the time, . Each bar element has a thickness (ℎ ), length ( ), inclination ( ) and a unit thickness in the out-of-plane direction.
The displacements over the element are assumed to vary linearly (linear shape functions) from which the corresponding strains (and therefore stresses) can be computed (more details in Appendix III). By choosing a small element size, the changes in inclination at the nodes are insignificant and therefore displacements and nodal forces are assumed continuous.

Constitutive equations
Both the sliding layer and weak layer constitutive relationships are described in =̄, ⋅ max wherē, is the average peak shear strength of element i of the sliding layer; is the engineering plastic shear strain required to soften the sliding layer material from peak to residual strength.
While other functions of strain softening may be appropriate, for example exponential, 45 Stöcklin et al. 31 have shown that the value of the shaded area in Figure 4 (and hence the choice of function) is of secondary importance to the postfailure geometry and that the value of residual strength plays a more dominant role.
Similar to the sliding layer, the undrained shear strength of the weak layer, , is defined by where , is the plastic displacement across the weak layer at node i, and is the plastic displacement required to soften the weak layer from its peak strength to its residual.

Analysis steps
Initial step-pre-failure stress state Because the elements used in the analysis have only a single degree of freedom, the pre-failure horizontal stresses cannot be obtained by simply ramping up gravity. Instead, the averaged horizontal stresses ( ′ ) are first assumed to correspond to infinite slope 0 stresses, that is where ′ is normal active stress given by Equation (35) in Appendix III. Gravity is ramped using a smooth step ramping and solved by an implicit solver, after which the initial stresses deviate from the 0 stresses until equilibrium is established.

Failure step
After the initial conditions have been calculated in the first step, the slide is triggered by an 'earthquake', represented by a pseudo-static horizontal acceleration applied by external nodal forces (see Appendix III for more details) for the duration of one second, after which the slide is driven by gravity alone. In this step, a dynamic explicit solution scheme is used to allow for faster computation time. In the elements experiencing significant stretching, an adaptive remeshing procedure is applied to reduce inaccuracies that arise from the elements' elongated dimensions. A remeshed element is halved, a new node representing the boundary between the two elements is introduced, and the element and nodal quantities are interpolated. Because the problem is highly nonlinear, the interpolation of nodal and element quantities may cause the explicit scheme to diverge since equilibrium is not ensured. For this reason, the calculation at this time increment is repeated using an implicit scheme to establish equilibrium before the explicit solution scheme is reused. This adaptive remeshing and solution schemes are the main reason why a finite element formulation was chosen over the more common finite differences method.

Solution update
The solution scheme of the dynamic problem is illustrated in Figure 5. The analysis time is broken into small time increments based on Equation (23), in which the nodal displacements and velocities are evaluated by where Δ is the nodal incremental displacements at time , −0.5 and +0.5 are the nodal velocities at times −0.5 and +0.5 , respectively, is the nodal acceleration while is the stable time increment with , the plane-strain modulus of element i of the sliding layer and is the sediments' total density. The new strain and stress increments are then calculated from the incremental displacements (see Appendix III). The nodal internal and external forces are assembled, and the nodal accelerations are calculated by solving the dynamic equilibrium equations of the system. Once all necessary quantities are computed and before calculating a new increment, the nodal positions are updated according to TA B L E 1 List of parameters for framework application.

Parameter Symbol Value
Half slope height (m) 40 Maximum inclination (degrees) 8.5 Thickness of the sliding layer (m) ℎ 10 Thickness of weak layer (m) 0.5 Plane-strain modulus (kPa) 3600 Shear modulus of weak layer (kPa) 1136 Total density (gr/cm 3  The new depth, , and inclination, are read based on the new nodal positions. The analysis is considered finished when the total kinetic energy of the system reaches a certain low threshold value ( Figure 5), that is, when all masses can be considered to have come to rest.

GIS BASED FRAMEWORK FOR IDENTIFICATION OF AREAS PRONE TO SUBAQUEOUS LANDSLIDES AND THEIR POST-FAILURE EXTENT
In the previous section, the theoretical background and methodology required to estimate landslide prone areas and their extent were described. In this section, the steps required to translate the input parameters (such as bathymetry and physical properties) and their implementation within the GIS framework will be elaborated. The suggested framework utilizes the aforementioned criteria and follows to some extent the process suggested by Rushton et al. 2 Our framework however incorporates few additional steps, which were not included in Rushton et al., 2 thus allowing to overcome some of the difficulties arising from their approach. In addition, the suggested framework provides the extent of upslope retrogressive failure and of confined failure downslope, which are critical for the risk assessment.
To ensure that the results are aligned with the theory presented earlier, the framework is first applied on an artificially constructed surface consisting of an extruded 2D exponential slope. For this sake, the parameters listed in Table 1 were considered, which closely follow those used by Zhang et al. 28 for their retrogressive failure study. After the framework is validated for this simplified case, it will be applied to the reconstructed bathymetry of the southeastern banks of the Küssnacht Basin in Lake Lucerne, central Switzerland, which failed during the 1601 earthquake. The complications arising from applying this approach to a complex 3D bathymetry will be discussed and possible solutions will be presented.

Dimensions of initial failure zones and the factor of safety
An extruded exponential slope with a total height-drop of 2 = 80 and maximum inclination of = 8.5 • was considered in the analysis ( Figure 6C). The initial slide dimensions can be assessed by determining whether the slope fulfills the failure criterion (Equation (10)). First, it is required to differentiate between the unstable, quasi-stable and stable zones based on shear stress ratio and measure the initial length of the unstable zone ( ). If the failure criterion is fulfilled, catastrophic failure is likely to take place and propagate to the external borders of the quasi-stable zones. The unstable zones in which the failure criterion is fulfilled are referred to as failure initiation zones, while the quasi-stable zones surrounding them are referred to as failure propagation zones.
In the following section, the steps required to derive the dimensions of the failure initiation and propagation zones as well as the FoS distribution using the GIS based framework are presented. Each step will be described and supported by complimentary figures illustrating the necessary process. 3.1.1 Unstable, quasi-stable, and stable zones The initial step involves calculating the spatial shear stress ratio distribution. This requires knowledge of the soil's physical properties and the spatial distribution of slope inclinations. The shear stress ratio distribution is calculated by applying Equation (6) to the obtained slope inclination distribution ( Figure 6B) together with the parameters in Table 1 using the built-in Raster Calculator tool. The resulting layer is a raster, containing cells with an attribute of the shear stress ratio value ( Figure 6D). The slope inclination raster ( Figure 6B) is derived from the bathymetry raster using the built-in Slope tool.
To distinguish between the unstable zone, where ≥ 1, the quasi-stable zone, where, 0 ≤ < 1 and the stable zone, where < 0, a reclassification of the shear stress ratio distribution is carried. Using the built-in Reclassify tool allows assigning a different integer value for each one of the zones ( Figure 6E). An integer raster is necessary, so that it could be converted into polygons designating the different zones. The polygon entities in ArcGIS are more flexible in terms of data processing and are used later to derive the distribution of FoS. After the integer stress ratio raster is generated, it is converted into polygons using the built-in Raster to Polygon tool.

Factor of safety
The shear band propagation approach determines whether a slope is prone to catastrophic failure based on the initial length of the unstable zone. If this length exceeds a certain critical length (Equation (10)) the slope will fail catastrophically and the FoS (Equation (1)) should be calculated with respect to the residual strength over the entire failure initiation and propagation zones. If, however, this length is not exceeded, the slope will remain stable, and the FoS should be calculated with respect to the peak strength everywhere. The initial length of the unstable zone ( ) must be, by definition, computed as the length of the line following the steepest path (i.e., normal to the contour lines) and extending to the boundaries of the unstable area. For this reason, the initial step ( Figure 7A) involves overlaying the extent of unstable and quasi-stable polygons ( Figure 6E) with the contour lines ( Figure 6A). The contour line representing roughly the mid-axis of the polygon (orange dashed line in Figure 7B,C) is seeded with equally spaced points. The slope gradient in each point is then calculated and a straight line is extended and cut by the original boundaries of the unstable polygon ( Figure 7C). This results in multiple lines representing the initial length of the unstable zone. The average length of the generated lines is then saved as an attribute. In the case of a simple extruded exponential geometry, all lines are of identical length. However, in a more complexed geometry, the shape of the polygons deviates from rectangles resulting in varying initial downslope lengths. In such case an average value is taken as will be shown later.
The critical length (Equation (10)) depends on the average shear stress ratiōand the mean inclination of the unstable zonē. Therefore, their values are calculated within the unstable polygon using the built-in Zonal Statistic tool and saved as an attribute. The initial length of the unstable zone is then compared to the critical length and the FoS distribution TA B L E 2 Resulting parameters comparison for factor of safety calculation.

Parameter
Analytical GIS within the unstable polygon and its corresponding quasi-stable polygon can be calculated according to Table 2 below summarizes the results of the calculations illustrated in Figure 7C and compares them to the analytical results. The calculations using the GIS framework show good agreement with the analytically derived values and will show better results with decreasing grid size. The calculated initial length is larger than the critical length and therefore the slope is susceptible to catastrophic failure due to the SBP. As a result, the factors of safety over the failure initiation zone and the enclosing failure propagation zone are both calculated using the residual strength. Figure 8 shows the obtained FoS distribution with the red lines illustrating the extent of the predicted failed areas. Figure 8A presents the obtained distribution of FoS according to the SBP approach. For comparison, Figure 8B shows the predicted FoS distribution based on LEM, which clearly underestimates the extent of the failure zones.

Growth of initial landslide dimensions due to tertiary failure mechanisms
So far, the methodology presented is capable of determining whether a given slope is prone to global failure due to catastrophic shear band propagation. The dimensions of subaqueous landslides can increase substantially outside the secondary failure zones due to tertiary failure mechanisms such as upslope retrogressive failure, ploughing, bulldozing, and runout. Therefore, the factor of safety distribution in Figure 8 is only indicative of the areas initially affected by catastrophic shear band propagation and global slope failure, and the final dimension of the landslide depends on the extent of the possible tertiary failure mechanism. Upon downward displacement of the failed slab, the intact soil at the head scarp position is unloaded. If the soil is unable to withstand the reduction in compressive axial load due to subduction of the failed slab, recurring retrogressive active failure takes place. Downslope the failed slab may either emerge and run out or alternatively remain confined and cause further shearing of the basin sediments below. In this section, the criterion for retrogressive failure of Zhang et al. 28 will be incorporated into the GIS framework and the procedure required to implement their methodology for estimating the retrogression length will be demonstrated. In addition, the DIFEM formulated earlier will be utilized to estimate both extents of the upslope retrogressive failure and of the downslope confined failure of the landslide. For this aim, the calculations will be carried out for the same extruded exponential geometry depicted in Figure 6A,C. Finally, the dimensions of different failure zones from all methods will be compared and discussed.

Retrogressive failure
According to Zhang et al., 28 evaluation of the extent of retrogressive failure requires calculation of the normalized initial sliding layer resistance (Equation (11)), the maximum normalized driving force (Equation (14)) and the stability number (Equation (13)). Table 3 below provides the calculated values. Since the normalized initial sliding layer resistance isˆ0 = 2.166 > 1, the case at hand can be classified according to Zhang et al. 28 as SBP dominant mechanism with overlaying layer of high strength. Because the sliding layer resistance is insufficient to withstand the maximum driving force (i.e., 1 <ˆ0 ≤ˆ), a slab failure will occur. With the stability number falling below the threshold, = 0.742 < 1, retrogressive failure will take place and its extent can be calculated using Equation (17). In the case of SBP dominant mechanism, the extent of retrogressive failure ( ) is defined by   28 as the length of spreading above the line where the current shear stress ratio is zero, that is, = 0 (see illustration in Figure 2), where is calculated using Equation (18). The calculation, however, requires an iterative process since the retrogression length in Equation (17) is a function of the average current shear stress,̄, calculated over the length .
The distribution is calculated and reclassified ( Figure 9C) using the slope raster derived from the slope geometry ( Figure 9A,B). The = 0 defines the line from which the retrogression length is found and measured. The distance of the = 0 from the middle of the slope is obtained at =0 = 358 m (see Figure 9C). To start the iterative calculation, an initial length of retrogression is assumed, , , after which the average current shear stress ratio over the assumed length, , is calculated. Then, Equation (17) is used to obtain the calculated retrogression length, , , usinḡ, . If the calculated length does not satisfy the criterion | , − , | ≤ , with being a certain small value threshold (chosen 1 m in this case), the solution has not converged, and a new length is assumed , ∶= ( , + , )∕2. The process is then repeated until the above condition is satisfied.
To illustrate the calculation process, an initial retrogression length of , = 100 was first assumed and the = 0 line was buffered by that amount (marked by a red hatched polygon in Figure 9D). Using the built in Zonal Statistics tool, the average value of within the buffered polygon,̄, , was calculated and found to be −0.0404. This value was then substituted into Equation (17) and the corresponding calculated length was derived , = 278.3 . A new assumed length was calculated by , ∶= ( , + , )∕2 = 189.5 and the process was then repeated by calculating the average value of within the new buffered polygon (marked by a blue hatched and striped polygon in Figure 9D) and calculating the corresponding calculated length. The process converged after four iterations where the obtained value of retrogression length was found to be 173.1 m. Hence, the final extent of the landslide in the upslope direction is found at =0 + = 531.1 from the mid-axis of the slope. Equation (17)  Solving the above equation yields a value of = 171.8 which is very close to that obtained using the suggested procedure implemented in GIS.

Downslope tertiary mechanisms
The DIFEM and its formulation outlined in Section 2.4 are a convenient and fast tool for determining the extent of subaqueous landslides. Both retrogressive failure and frontally confined failure can be modeled using this method. Theoretically, the method allows simplified modeling also of frontally emergent slides by releasing the stresses at the landslide front and updating the element dimensions accordingly. However, the problem of an emergent slide is significantly more complicated because of potential water entrapment, hydroplaning, turbidity currents etc. As a result, at this stage it has been decided to leave the problem of frontally emergent slides outside the scope of this paper. Extension of the developed DIFEM to account for frontally emergent slides is discussed later. Before applying the suggested numerical framework, a short study was conducted to validate its applicability. To this end, the commercial finite element software ABAQUS 47 was employed. The coupled Eulerian-Lagrangian (CEL) technique, built within ABAQUS computing environment, allows materials to undergo very large deformations due to the Eulerian implementation. This way, the problems related to mesh distortion associated with standard Lagrangian finite element analysis are minimized. The CEL method was applied in various studies to model similar problems and in back analysis of past slides, for example, Dey et al., 30,48 Stöcklin et al., 31 Klein and Puzrin, 19 Zhang et al., 28 Klein et al. 32 Similar to DIFEM, the slide is triggered using a pseudo-static horizontal acceleration applied over 1 s, after which the slide becomes self-driven by gravity only. More elaborate details of the modeling technique can be found in Trapper et al., 49 Stöcklin et al., 31 Klein and Puzrin 19 and Klein et al. 32 Figure 10 shows the post-failure slide geometry and the distribution of equivalent plastic strains as obtained by the ABAQUS model. The dashed orange line represents the post-failure geometry as obtained by DIFEM presented in this paper. In addition, the extent of retrogression calculated earlier is also shown. Upslope, Zhang et al. 28 methodology and DIFEM provide similar and slightly more conservative results (about 30% larger) in terms of the retrogressive failure extent compared to the ABAQUS model. Downslope, DIFEM also shows 15% larger extent ( in Figure 10) of displaced material in comparison to the ABAQUS model. The difference originates most likely from the depth-integrated method consisting of single sliding layer elements in which internal shear planes cannot develop like in the ABAQUS model (see the zoomed in frame in dashed red line in Figure 10). As a result, less energy is dissipated in shearing of the sliding material, leading to increased kinetic energy in DIFEM compared to the ABAQUS model. The increased kinetic energy results in further displaced material, which ultimately may lead to conservative estimation of the landslide extent. From this comparison it can be concluded that both Zhang et al. 28 methodology and DIFEM are useful and adequate tools for predicting conservative ultimate landslide dimensions.
The results obtained from the above analyses were projected onto the map illustrating the susceptible area depicted in Figure 8A and the resulting dimensions of the tertiary failure are shown in Figure 11. As expected, the downslope boundary of the landslide extends beyond the initial landslide dimensions predicted by the Puzrin et al. 15 criterion. It is interesting, however, that both FE methods applied, and the procedure suggested by Zhang et al. 28 resulted in upslope failure boundary which is confined within the area prone to catastrophic failure suggested by Puzrin et al. 15 This is likely due to the difference in which the earthquake induced shear stress is considered. Puzrin et al. 15 shear stress ratio considers a constant shear stress applied by the earthquake loading. Depending on the magnitude of the horizontal acceleration, even relatively flat zones may be considered quasi-stable, especially when the earthquake induced shear stress is considerable comparing to the gravitational shear stress. Although relatively flat zones may suffer some displacements during an earthquake event, once the earthquake ceases, there may be no additional sufficient driving force that will cause any further displacement and hence catastrophic failure involving these regions. The retrogressive failure methodology of Zhang et al. 28 evaluates the retrogression length after the triggering process has taken place. Similarly, both FEM analyses introduced earlier apply the horizontal acceleration for a duration of 1 s only, after which the slide moves solely due to gravity. As a result, the post-trigger quasi-stable zone is much smaller than that obtained by Puzrin et al. 15 This suggests that the criterion of Puzrin et al. 15 may produce conservative estimates of the affected areas upslope. However, since more studies are required to investigate the possible extent of retrogressive failure upon seismic shaking, it is impossible to say in this moment that applying the Puzrin et al. 15 criterion will always provide a safe estimate. As a result, it is suggested to still consider either Zhang et al. 28 or DIFEM to estimate the extent of retrogressive failure.

APPLICATION TO LAKE LUCERNE AND THE ZINNEN SLIDE
The suggested framework has so far been applied to an artificially generated slope surface. In this section, we demonstrate its application to the reconstructed bathymetry of the Küssnacht Basin, Lake Lucerne, Switzerland. The goal of this application example is to illustrate some of the difficulties that arise by applying criteria that were formulated for 2D slopes onto 3D bathymetry and to present some possible solutions. Moreover, the application to the well-documented Zinnen slide and its surrounding will demonstrate the framework's ability to predict subaqueous slope failures and their final size.

Study area-southeastern banks of the Küssnacht basin
The Küssnacht basin is one of seven sub-basins of Lake Lucerne. It is characterized by slopes reaching 35 o and nearly flat deep-water basin-plains. 1,50,51 The sedimentological sequence consists of fine-grained Holocene sediments overlying Late Glacial and Glacial deposits. The Holocene sediments are about 4−5 m thick on the slopes and thicken to about 5−15 m at the lake basin floor. 52 During an earthquake around the year 1601, a significant portion of the Küssnacht basin slopes failed 52,53 causing eventually a tsunami wave of about 4 m that flooded the city of Lucerne. 54 The slides were shown to take place over the Late Glacial and Glacial deposits interface. 53,55 Recent studies investigated the failure of the Zinnen Slide on the southeastern banks of the Küssnacht Basin ( Figure 12A). A field campaign was carried to characterize the Zinnen Slide where a high-resolution seismic reflection profile was obtained, free-fall CPTU tests were performed, and undisturbed samples of the upper Holocene sediments were retrieved and tested in the lab. 53 The Zinnen Slide area shows thinner and stronger sediments on the slopes compared to those in the basin ( Figure 12B,D), which contributed the Zinnen Slide to fail in a frontally confined manner. 32 The sediments' strength increases linearly with depth and the strength of the weak layer was found to be a fraction of the strength (denoted as ) of the overlying sliding sediments. Table 4 lists additional parameters required for the current study (more details in Klein et al. 32 ).
In this study, the bathymetry of the Küssnacht Basin's southeastern banks was reconstructed. The contour lines of the existing bathymetry were modified in the areas of failed slopes based on constant sediment thickness assumption. 57 Using the new contour lines, a new digital elevation model was created on which the framework was applied.

Model calibration
Establishing the unstable and quasi-stable zones, based on the shear stress ratio (Equation (6)), requires determining the site-specific pseudo-static acceleration coefficient multiplier ( -Equation (4)). Dimmock et al. 58 describe a procedure for determining this coefficient for the West Nile Delta region based on Newmark analysis. By permutating the input parameters, they determined the coefficient corresponding to a factor of safety of unity (based on LEM calculation) for varying permanent Newmark predicted displacements. By setting a permanent displacement threshold of 0.1 m, they determined the appropriate average pseudo-static coefficient multiplier of = 0.23 with a standard deviation of 0.08.
Here we suggest a different site-specific approach for calculating the coefficient. For this aim, we seek a critical ℎ value that brings a reconstructed failed slope to failure based on the SBP approach. Furthermore, the initial dimensions of the failed mass should match well those observed. The value of ℎ was incrementally increased from zero such that in each Thickness of sediments in the basin (m) ℎ 9 Stability number *

16.9
Shear strength increase constant of weak layer * 0.413 *see Table II.1 in Appendix II for derivation.
iteration, the length of the unstable zone increases (based on Equation (6)) until it meets the critical length (Equation (10)) in which case the slope was considered to have failed due to SBP. After the critical ℎ is obtained, the value of is back calculated using Equation (4).
To establish the appropriate value for the Küssnacht Basin, the suggested approach is applied to the Zinnen Slide. For this aim, the pre-failure geometry of the slide was reconstructed, and the parameters listed in Table 4 were used. Figure 13A,B respectively, show the current and reconstructed hill shade and contour lines of the Zinnen slide area.
The process begins with Figure 13C showing the derived unstable and quasi-stable zones for the case of ℎ = 0, that is, without any pseudo-static seismic loading. In this case, the model must show that the slope is stable. The model indeed results in initial unstable zone length of = 86.5 m, which is smaller than the calculated critical length of = 146.5 m, reassuring that the slope is stable under static conditions. For calculating the critical length (Equation (10)) test data concerning the ratio ∕ is required. In absence of this data, a conservative assumption is considering ∕ = 1, which was also used in the numerical analysis of the slide event. The obtained critical length was calculated using Equation (10), wherēand̄were derived using zonal statistics of the unstable zone. Figure 13D shows the derived failure initiation and propagation zones for the critical value of ℎ ≈ 0.17 under which the slope is considered to have failed. For this value of ℎ the average initial length was found = 116.2 which exceeds the calculated critical length of = 92.1 . The area covered by the failure initiation and propagation zones agrees well with size of the failed portion of the Zinnen slide (excluding the extent of tertiary failure mechanisms, discussed below). With peak ground acceleration of 0.14 g, 57 the resulting pseudo-static acceleration coefficient multiplier for the Küssnacht Basin is approximately = 0.263.
Using the derived value of ℎ the framework was applied to additional parts of the Küssnacht Basin where the pre-failure bathymetry was reconstructed. Figure 14 shows the final results of the analysis conducted on the southeastern banks of the Küssnacht Basin. The SBP approach used here to estimate the safety of subaqueous slopes was derived for slopes that can be described by plane-strain conditions. Under these conditions, the shear band grows along the steepest gradient line, that is, normal to the contour lines. To accommodate these limitations, the slope is divided based on its aspect into segments representing near plane-strain conditions. To this end, an aspect analysis is conducted on the bathymetry raster. Aspect analysis assigns each pixel of the bathymetry raster with a value indicating the compass direction of the slope face. The values vary from 0 to 360 such that 0 represents north and the aspect value increases clockwise. Because the slope geometry is complex, the resulting aspect raster is reclassified using a certain reclassification interval. This interval is a user-defined variable that needs to be chosen based on the complexity and general orientation of the analyzed bathymetry. In this case a threshold of 45 • was chosen. After the aspect raster is obtained and reclassified, it is converted into a polygon feature class using the built-in Raster to Polygon conversion tool. This results in polygons representing the different slope aspects 1−11 in Figure 14. The failure initiation polygons (2, 3 and 11) fulfill the failure criterion, that is, ≥ . Polygon 5 (marked in thick dashed line) is a border case, since it shows an initial length, , lacking only 6% its value to exceed the calculated critical length, (87.1 m as opposed to 92.4 m). The rest of the polygons did not fulfill the failure criterion by a large margin (values in Table II.2 in Appendix II) and are therefore considered stable. The polygons 6−10 which were found stable agree with the current bathymetry, that is, they indeed did not fail during the 1601 earthquake. The only polygons that do not agree with the observed geomorphology are polygon 4 and partially polygon 1. They represent, however, less than 5% of the total analyzed area. The discrepancy may be due to the reconstruction procedure of the pre-failure geometry and the aspect-based splitting of the polygons. Their possible impact is discussed in the limitations section below.

Final slide configuration
The Zinnen Slide pre-failure geometry ( Figure 15) was derived by creating a cross-section along the seismic line in Figure 12C (section A in Figure 14). The thicknesses of the sediments on the slope and in the basin were interpreted based on seismic reflection profiles by Sammartini et al. 53 (more details in Klein et al. 32 ). Figure 15 shows the variation in strength profiles of the slope and basin sediments derived based on CPT results and FEM analysis. 32 The geometry and strength input parameters were used to run DIFEM calculation of the Zinnen Slide and the results, in terms of postslide geometry, were compared with CEL-FEM analysis and the interpretation of seismic reflection profiles. 32,53 Figure 16 F I G U R E 1 4 Failure initiation (primary failure) zone ( ≥ 1) and propagation (secondary failure) zone (0 ≤ < 1) of the southeastern banks of the Küssnacht Basin. Sections calculated in DIFEM (thin black lines) and the extent of the calculated slides (thick yellow lines).

F I G U R E 1 5
Pre-failure geometry of the Zinnen Slide after Klein et al. 32 presents the results of the slide's post-failure geometry where both CEL-FEM and the proposed DIFEM ( Figure 16A,B) show very good agreement with the observed geomorphology on site ( Figure 16C), with DIFEM providing a reasonably conservative estimate (about 10% larger measured from the toe). Aside for the Zinnen slide, two more cross sections were calculated (section B and C in Figure 14) using the DIFEM. The sections geometry was smoothened to allow the DIFEM calculations. Once the results were obtained, the upslope and downslope extents were extracted and projected onto the cross-section lines to illustrate the predicted area affected by the slide marked in thick yellow lines in Figure 14. When compared to the corresponding border of the thrust and fold belts (white solid lines in Figure 14 as interpreted by Sammartini et al. 53 ) all three calculations provided good fit to the observed geomorphology. The calculation results for cross-sections B and C are provided in Figure II.1 in Appendix II. The southeastern banks of the Küssnacht Basin do not show any physical evidence of retrogressive failure. The analytical criterion for retrogression, that is, Equation (15), has also not been fulfilled (see Table 4). Therefore, the methodology for calculating the retrogression length was not employed. Both the DIFEM and the CEL-FEM calculated using ABAQUS agree on the upslope extent of the slide and support the analytical criterion suggesting that no retrogressive failure could be expected ( Figure 16A).

LIMITATIONS, CHALLENGES, AND FURTHER RESEARCH
Submarine and sub-lacustrine environments consist of complex geomorphology and sedimentary sequences. Consequently, the application of criteria and frameworks that were formulated for plane-strain conditions and constant sediment properties and geometry raises difficulties in analyzing the dataset. As a result, some simplifications and assumptions are required.
In this framework we applied the criterion of Puzrin et al. 15 as it is simple to implement, conservative and considers 2D geometry. As a meaning to apply it to 3D geometry, the use of aspect analysis was required. The choice of the aspect threshold and simplifications of the polygons' borders should always be overlaid with the contour lines to ensure meaningful division. Moreover, the way the polygons are divided influences the calculated values of the initial and critical lengths of the unstable zones and consequently the determination of failure. For example, polygon 3 in Figure 14, when divided differently, results in a smaller failed portion (see Appendix II and Figure II.2 for more details). Therefore, it is recommended to attempt dividing the polygons in such a way that they have roughly the same initial length. This way, the influence of using a mean length value of equally spaced lines (see Figure 7) is also minimized.
Conforming with the SBP approach, we assume that failure takes place along the steepest gradient lines. This assumption may have its limitation when considering a margin polygon (e.g., see polygon 1 in Figure 14). In this work, we use the outermost lateral extent of such polygons and construct lines normal to the contour lines to create the associated quasi-stable polygons, 18,19,42 however, have studied the effect of out-of-plane geometry on the critical shear stress ratio for planar and conical slopes. Their studies showed that the shear surface may grow in the transverse direction due to downslope movement of the sediments and that the growth extent is related to the aspect ratio and area size of the pre-softened zone (or failure initiation zones in this case). The proposed methodology does not take this effect into consideration and therefore may result in non-conservative estimates around these specific areas. Therefore, more effort is encouraged to investigate the implementation of failure criteria that extend beyond the plane-strain assumption. Additional geomorphological shapes that require further study are berm-like slopes, domes, and valleys.
For back-analysis of past slides, reconstruction of the pre-failure geometry is required, which introduces uncertainties. Therefore, the results should be adopted with great care and an engineering judgment must be used as the slope steepness has significant influence on the failure estimation. That said, forward-analysis using the suggested framework to existing and intact slopes does not suffer from this limitation.
In the analysis of the Küssnacht Basin, we used uniformly distributed physical properties to describe the sediments' behavior based on tests conducted at the Zinnen Slide site, although they might vary across the basin. The suggested framework, given a better knowledge of spatial distribution of physical properties, can easily account for such variations by creating raster datasets and using them instead of constants in all the necessary formulas.
The digital elevation model raster cell size and respective contour intervals should also be considered. If the raster cell size is comparable or significant compared to the measured and calculated lengths, the results should be handled carefully.
The developed DIFEM presented in this study currently can only be applied for confined slides. Presently, it cannot accommodate the change in material behavior once the sliding sediments can no longer be described by an incompressible solid as in the case of frontally emergent slides and turbidity currents. In cases where mixing with water is limited, the suggested DIFEM can be extended to model mudflows by accounting for a more appropriate constitutive behavior such as Bingham fluid.

SUMMARY AND CONCLUSIONS
In this paper we presented a comprehensive GIS framework aimed at analyzing deterministically the stability of submarine and sub-lacustrine slopes. The proposed methodology considers an advanced failure criterion, based on the SBP approach, 15 which accounts for the strain-softening nature of soft clay-rich marine and lacustrine sediments. Application of this criterion allows determining the initial dimensions of the susceptible unstable areas using the physical properties of the sediments, available bathymetry in the form of digital elevation model, and a triggering process (a pseudo-static earthquake load). After the primary and secondary failure areas are determined, the suggested framework provides tools to estimate the final, post-failure configuration of the slide both at the crown and at the toe. An iterative process for calculating the distance with which failure retrogresses upslope, based on the principles of Zhang et al., 28 was presented. A mixed implicit-explicit dynamic depth-integrated finite element method (DIFEM) with adaptive remeshing was developed to estimate the post-slide geometry (in the case of frontally confined slides). Both iterative retrogressive failure process and DIFEM were benchmarked against large deformation finite element analysis (ABAQUS CEL-FEM), showing good agreement or slightly conservative results. After benchmarking the suggested framework against a simple case of an extruded exponential slope, it was applied to the documented case of the Zinnen Slide at the southeastern banks of the Küssnacht Basin in Lake Lucerne, central Switzerland. For this aim, the pre-failure geometry of the failed slopes was reconstructed. A procedure for choosing the appropriate pseudo-static earthquake loading coefficient based on the SBP approach and an emphasis on agreement with the slide's dimensions was illustrated. Then, the framework was applied to the Küssnacht Basin and showed promising results in predicting the slope failures that took place during the 1601 earthquake. The developed DIFEM was then applied to the Zinnen Slide, where a seismic reflection profile is available, and showed good agreement with the observed geomorphology of the slide and with CEL-FEM simulations. Additional DIFEM calculations across the Küssnacht Basin show good final landslide dimensions prediction.
The use of analytical solutions and simplified numerical methods (e.g., DIFEM) can be questioned when commercial codes such as ABAQUS's CEL-FEM with their great capabilities are available. This is indeed a legitimate question for calculation of an individual slope. For the probabilistic risk assessment of entire basins, however, it is not feasible to run millions of such calculations due to their long runtime and computationally expensive resources demand. An attempt should be made to develop reliable analytical solutions, which can be easily incorporated into a GIS-based framework. If such analytical solutions are unavailable, the effort should be invested into development of simple and computationally inexpensive numerical frameworks (e.g., depth-integrated methods).
While the proposed framework provides an automated procedure to overcome the need for the manual intervention reported by Rushton et al., 2 the results should still be reviewed in light of the limitations of the framework as discussed in the limitations section and in Appendix II. Guidelines for dealing with these limitations are also provided in the paper.
Within this context, the main findings of this paper can be summarized as follows: • The shear band propagation approach appears adequate for analyzing the stability of subaqueous landslides. Its application and incorporation into GIS framework require additional steps (detailed above) comparing to simple limit equilibrium calculations, however, the limit equilibrium method falls short when applied to strain-softening sediments. • Analytical calculation of the retrogression length using the proposed methodology based on the work of Zhang et al. 28 provides good estimates when compared to simulations by large deformation finite element analysis. Attempts should be made to validate it against real-life cases in the future. • The depth-integrated finite element method is a useful, fast, and suitable tool for calculating landslides' post-failure geometry in the case of confined slides. Its formulation should be extended to allow simple calculations of emergent runout distances.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The data that supports the findings of this study are available in the manuscript and its appendices. Supporting data concerning the analysis of the Küssnacht Basin can be found in Sammartini et al. 53 3. The starting point of retrogression can be found by setting the current shear stress ratio to zero and solving for the coordinate : 4. The retrogression length is given by Equation (17) in which̄is the average current shear stress ratio over the length Which can be calculated analytically by The derived value is inserted into Equation (17), which after simplifications results in The above equation can be solved numerically to obtain the retrogression length.

APPENDIX II: ANALYSIS OF THE KÜSSNACHT BASIN
Results of the additional cross-sections B and C using DIFEM Figure II.1 presents the results of the DIFEM analysis of cross-section B and C from the Küssnacht Basin (see Figure 14). In both models, the geometry was extended by 50 m upslope (negative coordinates areas) to minimize boundary effects. Downslope, the models were extended even further, however their full extent is not depicted in the figure for illustration purposes. Due to lack of seismic reflection profiles, it was assumed that as in the Zinnen slide section, the sediment thickness on the slope prior to the collapse was 6 and 9 m in the basin. (Tables II.1 and II.2)

The influence of aspect polygon division on determination of failed slopes
To illustrate the possible effect of aspect polygon division on determining whether a slope is prone to failure, Figure 14 illustrating the derived failure initiation and propagation polygons was redrawn using a different unstable polygons division ( Figure II.2 below). Specifically, polygon 3, which shows large variation in the initial length of the unstable zone, was divided into two separate polygons 3a and 3b, to create polygons with more evenly measured initial length. When recalculating both initial and critical lengths of the two new polygons it appears that only polygon 3a fulfills the failure criterion (Equation (10)) with = 116.6 < = 108.2 while polygon 3b becomes stable with = 67.8 < = 101.0 . Consequently, the aspect polygon division is a useful tool to divide the complex geometry into 2D plane-strain like polygons. However, the division must be inspected in order to avoid such cases as discussed above.

APPENDIX III: DIFEM FORMULATION
After calculating the incremental displacements, the strain increment in the sliding layer is defined as The normal, slope parallel stress ( ′ -positive in compression) and strain ( -positive in compression) are calculated at the integration points (see Figure 3). An elastic predictor plastic corrector method (Figure 4) is employed to calculate the stresses in the sliding body using the obtained strain increments. It is first assumed that the strains are fully elastic, and the new stress increment is calculated The average slope-parallel stresses in the sliding layer are expressed in terms of the average undrained shear strength bȳ′ wherē′ , ,̄′ , ,̄′ , are the average normal effective stress, active failure stress and passive failure stress, respectively and , ( , ) is given by Equation (19). The new stresses are then compared with the current active and passive failure stresses ′ −1 , Wherē′ −1 , ,̄′ −1 , are calculated using Equation (35). Because the sliding layer is modeled using bar elements where no shear strains can be computed, it is assumed that the stresses and strains in the sliding layer are aligned with the principal directions and therefore Where, Δ , and Δ , are the maximum and minimum principal strain increments and Δ , is the strain increment in the normal direction. The incompressibility condition (under plane-strain assumption) is expressed by Which therefore together with Equation (37) result in It is assumed that no slip occurs between the weak layer and the stiff base below it. Therefore, the shear strain increment in the weak layer is defined by where is the weak layer thickness. The current stress in the weak layer is updated using the new nodal displacement and a similar procedure of elastic predictor plastic corrector as in the sliding layer but using the corresponding constitutive law (Figure 4 and Equation (20)). Based on the incremental displacement, the new element dimensions are calculated using where, 0 is the element's initial length and is the element volume. Note that the element volume generally does not change throughout the course of the analysis unless an element is remeshed.
The nodal internal forces, , , are in equilibrium with the element stress, that is The weak layer is considered in the analysis as an external force acting at the nodes. The mobilized shear stress in the weak layer, , , is calculated in each node based on its respective displacement and the constitutive relationship depicted in Figure 4 and Equation (20). The external force applied at the node, , , is then calculated as To calculate the external gravity load, , , the mass of the element is distributed equally to its nodes such that where is the water density. To trigger the slide, a pseudo-static earthquake load is applied in the horizontal direction. Similar to the gravitational body force, the earthquake load, , , is distributed equally to the element nodes as follows where ℎ is the applied horizontal acceleration. After the stresses in both sliding layer and weak layer are obtained, the internal and external nodal forces are assembled, and the residuals, , are calculated by Where , is given by Equation (42) where and are the nodal acceleration and mass.
In every time step, the kinetic energy of the system is computed by Equation (48). If the kinetic energy reduces and drops below a certain threshold ( Figure 5) the analysis is considered completed.
where is the number of nodes in the analysis.