Kinematic energy balance approach to submarine landslide evolution

The paper applies the energy balance kinematic method of plasticity theory to the large deformation problem of initiation and propagation of the spreading and ploughing failure outside a failed slab in submarine sediments. The models account for the phenomenon of the progressive propagation of a slope parallel slip surface, which is also quantified using the energy balance approach. In contrast to existing approximate analytical and numerical solutions, the proposed approach provides a theoretical basis for spreading and ploughing criteria as well as the comprehensive dynamic solution of the problem of post-failure landslide evolution. Incremental integration of the derived analytical expressions for kinetic energy in time allows for modelling recurrent initiation of new kinematic failure mechanisms with their subsequent large-scale deformation. Treating the failed slab as well as the spreading and ploughing mechanisms as one composite dynamically evolving mass movement allows for the final post-failure geomorphology of the failed slope to be predicted using basic mechanical principles.


Introduction
Gigantic dimensions and peculiar geomorphology of historic submarine landslides have generated significant scientific interest in the literature (e.g., Hampton et al. 1996;Locat et al. 2002;Masson et al. 2006), in particular due to their potential tsunami hazard (see recent reviews by Grilli et al. 2017;Lovholt et al. 2017). Offshore and coastal developments have further highlighted the need for reliable hazard assessment, requiring understanding the mechanisms of the dynamic landslide evolution -from pre-failure conditioning, triggering, nucleation and propagation of localized shear zones to global failure and post-failure processes.
The shear band propagation (SBP) approach for submarine slope stability analysis helps to explain how a relatively short initial slip surface (localized shear zone of finite thickness, marked in red in Figure 1) can, under specific conditions, propagate dynamically tens of kilometres parallel to the slope surface through the quasi-stable part of the slope (orange in Fig. 1), eventually causing a slab failure. Proposed by Puzrin and Germanovich (2005) based on Fracture Mechanics principles developed by Palmer and Rice (1973), this approach allows determining the landslide geometry (for a slab failure), volume and breakoff velocity, providing more reliable input parameters for debris flow and tsunami wave propagation analysis (Puzrin et al. 2010;Trapper et al. 2015). In the past few years a number of novel analytical, numerical and experimental tools have been developed and validated (Viesca and Rice 2012; Puzrin et al. 2015Puzrin et al. , 2016Germanovich et al. 2016;Zhang et al. 2016Zhang et al. , 2017, extending the SBP approach to account for actual non-linear slope geometries, different types of triggers (excess pore water pressures, seismic loading) and dynamic propagation effects. Moreover, an attempt has been made to relate the observed spreading and ploughing mechanisms (formed above and below the failed slab, respectively, Fig. 1) to progressive propagation of a slope parallel slip surface into the stable part of the slope using simplified analytical criteria (Puzrin 2016;Puzrin et al. 2017).
D r a f t 4 In parallel, a significant body of work has been devoted to studying spreading and ploughing problems numerically. Progressive spreading failure accompanied by an uphill growth of the shear band in sub-aerial landslides driven by the removal of the downslope support has been studied using FE analysis (e.g., Locat et al. 2011;Locat et al. 2013;Quinn et al. 2011;Quinn et al. 2012;Dey et al. 2015). In more recent publications (e.g., Trapper et al. 2015;Dey et al. 2016;Zang et al. 2018) large deformation FE analysis was applied to explore the SBP driven post-failure evolution of submarine landslides. While sophisticated finite element analysis can provide insight into the basic landslide mechanisms its large computational costs limit hazard assessment to deterministic analysis of small seafloor areas (e.g., Locat et al. 2002;Vanneste et al. 2013).
Probabilistic slope stability analysis of large marine basins requires therefore an analytical framework not only providing criteria for spreading/ploughing initiation, but also predicting the post-failure deformations. In this respect, particularly promising is the energy balance approach proposed by Kvalstad et al. (2005) and used by Lovholt et al. (2016) for modelling the kinematics of submarine spreading failures. Studying gravity tectonics on passive margins, Peel (2014) suggested that the energy approach should be applied to the entire mass movement, including the failed slab, spreading and ploughing regions. The present work combines the energy balance kinematic approach of plasticity theory with the energy balance approach to the SBP, providing a more rigorous analytical framework describing the large deformation landslide phenomena as a synthesis of various interconnected post-failure processes.

Problem formulation
The phenomena explored in this paper take place after a slab failure of a portion of submarine slope. The movement of the slab is preceded by the formation and propagation of a slope parallel slip surface (SPSS) at its base. Due to this SPSS propagation, the intact slab experiences loading in its downslope and unloading in the upslope portions, respectively, which can eventually reach compression and extension failure within the slab body. Once failed, the slab often runs out at the D r a f t 5 bottom and leaves a scarp at the top of the slope. Under certain conditions, however, the slope continues to fail both up-and downhill from the intact slab (Fig. 2), with the kinematic failure mechanisms recurring within the sliding body. This paper explores these conditions and the subsequent failure evolution using rigorous energy based kinematic methods. The solution of this problem is achieved through the following three steps: • Formulate conditions for the initial slab failure; • Describe the evolution of the first failure mechanism in the sliding body, accompanied by large deformation and softening; • Formulate conditions for the recurrence of the failure mechanisms up-and downslope (spreading and ploughing).
The analysis of these three stages is performed with the help of a free body cut separating the intact slab from the rest of the slope (Fig. 2). This allows for the failure evolution at the two slab edges to be investigated independently. Interaction of the intact slab with the failed edges takes place through the slope parallel forces ܳ ௨ and ܳ , at the top and at the bottom, respectively. The criteria will be developed to determine when these forces cause initial failure and at what level they can cause subsequent retrogressive spreading at the top and progressive ploughing at the bottom of the slope. While the original intact slab is curved, the regions where failure mechanisms develop are assumed (for simplicity) to have a constant slope inclination.
After the up-and downslope failures are described separately, it will be demonstrated how the approach can be applied to the whole curved slope as one mechanism, allowing for defining the final geomorphology of landslide deposits.

Assumptions
The solution of the above problem requires certain assumptions with respect to the material behaviour. Due to the low permeability of submerged marine sediments and the fast loading rates, purely undrained failure conditions are assumed. The contractive behaviour of loose marine D r a f t 6 sediments often leads to strain softening after reaching a peak shear resistance. In this study the complex undrained soil behaviour of the sliding body is idealized by assuming linear elasticplastic behaviour with the Tresca failure criterion, where the undrained shear strength ‫ݏ‬ ௨ ሺߛ ሻ softens with increasing plastic shear strain ߛ from its peak ‫ݏ‬ ௨ to residual value ‫ݏ‬ ௨ (Fig. 3a). In case of strain localization within the sliding layer (formation of a shear band), the displacement jump ߜ along such a shear band is related to the plastic strain and the width of the shear band ݀ ௌ by the expression Although in reality the undrained shear strength varies in depth as a function of the initial effective stress and degree of overconsolidation, here, as a first approximation, a uniform strength distribution (averaged over the thickness of the sliding layer) is assumed. It is also assumed that the sliding layer is bounded at the bottom (at the constant depth h) by a weaker layer with different properties. Recent site investigations have indicated that submarine landslides are typically initiated at a weakened layer, within which the ratio of undrained shear strength to the effective normal stress is lower than for adjacent layers Gee et al. 2005;Sultan et al. 2010;L'Heureux et al. 2014;Gray et al. 2015). The elastic response of the sliding layer in slope parallel direction is important for the SBP along the SPSS. Following Palmer and Rice (1973) and Puzrin and Germanovich (2005), this elastic deformation is taken into account for the assessment of the propagation length of the shear band. In contrast, the threshold force at which the slab starts to develop plastic deformation (that is, forms a mechanism) at its edges and the subsequent evolution of the plastic deformation are formulated with the assumption of rigid-plastic strainsoftening material behaviour. This assumption is justified by the fact that in the failed region of the slope outside of the slab, the neglected elastic strain energy rate is small compared to the rate of plastic dissipation and the change in the potential energy.
D r a f t 7 Following (Puzrin and Germanovich 2005), the soil within the SPSS is assumed to be rigid-plastic with the peak strength value ߬ and exhibits strain-softening behaviour until the residual value ߬ is reached at the characteristic displacement ߜ (Fig. 3b). The deformations below the SPSS are neglected.
Although loose marine sediments usually show an increase of density with effective stress (close to the ground surface), for the sake of simplicity the total unit weight of soil ߛ is assumed constant in this paper. For submerged landslides only the effective unit weight of the soil ߛ ᇱ = ߛ − ߛ ௪ contributes to driving forces. In contrast, inertia forces depend on the total soil density ߩ = ߛ ݃ ⁄ .

The kinematic energy balance (KEB) approach for large deformations
After the initially intact slab body fails in extension at the top and compression at the bottom of the slope, it keeps unloading/loading the adjacent slope regions by the drop/heave in the seabed surface due to large deformations of the respective failure wedges. During this process the geometry of the failure zones changes dramatically and therefore the effects of not only material but also geometric softening and hardening have to be taken into account.
Accounting for the complete phenomena of elastic-plastic loading or unloading with geometric as well as material hardening and softening proves to be a complicated task, that in general can only be solved numerically using, for example, the finite element method. In order to develop a simplified analytical approach to the analysis of post-failure evolution in submarine landslides an extended kinematic method of limit analysis can be applied to describe plastic deformations (e.g., Kvalstad et al. 2005). The method has certain similarities to the upper bound-based strain path method proposed by Einav and Randolph (2004) and to the method applied by Alonso and Pinyol (2010) for the evolution of the Vaiont landslide. Similar to the classical kinematic method (e.g., Drucker and Prager 1952;Chen 1975) a kinematically admissible velocity field (i.e., a collapse mechanism) is assumed, based on which the work equation is formulated. In contrast to the classical kinematic method, however, where the work equation is formulated based on a virtual D r a f t 8 velocity field, here finite velocities and displacements are considered. Consequently, in order to fulfil the kinematic conditions (i.e., mass conservation) and boundary conditions, material has to be sheared across velocity discontinuities leading to a change in geometry. Similar to the classical kinematic method, the most critical collapse mechanism can be found by optimization of the geometry. Strictly speaking this should be done for every deformation state of the failure evolution.
In the present study, however, it is assumed that the optimized location of the velocity discontinuities for the initial collapse mechanism (i.e., corresponding to zero displacements) is also valid for the subsequent large failure deformation. Indeed, once the initial strain localization took place, due to the material softening it is highly likely that shear deformations will continue to be confined to the existing velocity discontinuity. The current state of the system is obtained by integrating the energy rates over time starting from the initial state. This approach allows not only taking into account the geometric nonlinearity, but also material softening, since finite strains and displacements can be calculated by integrating velocities. Furthermore, if the system is unstable, this approach allows describing the dynamic post-failure evolution by calculating the total kinetic energy of the system at every state of the loading history. In the following sections this method, which in this paper is denoted as kinematic energy balance (KEB) approach, is applied for assessing the post-failure evolution of submarine landslides, first at the top (unloading) and then at the bottom (loading) of the sliding slab.

Initial failure due to unloading at the top of the sliding slab
In this section the KEB approach is applied with the goal of formulating the condition for initial failure at the top of the sliding slab (onset of the development of plastic deformation in the slab body) and the subsequent evolution of the force in the upper part of the slope as a function of the state of plastic deformation. In order to achieve these goals, first, a kinematically admissible mechanism is introduced and the rate of kinetic energy is calculated as a function of the state of deformation. Then, following the concepts of limit analysis, the geometry of the mechanism for the D r a f t 9 initial failure state is optimized using the principle of minimal energy. Finally, the condition for the failure initiation is formulated in terms of the limiting slope parallel force.

Initial unloading failure mechanism
The proposed kinematic mechanism (Fig. 4) has one free geometric parameter ߱. It consists of three blocks and two velocity discontinuities VD1 and VD2 inclined, respectively, by the angles ߱ and ߱ + ߨ 2 ⁄ with respect to the slip surface. Block 1 represents a part of the sliding slab of thickness ℎ, which moves downslope with velocity ‫ݒ‬ ଵ . The remaining part of the slab is replaced at a free cut x by its slope parallel reaction ܳ ௨ , which will contribute to the work of external forces together with the gravity forces. Block 2 is an active wedge moving with the velocity ‫ݒ‬ ଶ parallel to the velocity discontinuity VD2, while Block 3 remains stable (conditions for the propagation of failure below and across the Block 3 are derived later in this paper).
The displacements of the deformed system are found by integrating the velocities of the respective blocks: The velocities and displacements of Blocks 1 and 2 are, therefore, both related by The corresponding drop in height of Block 2 is In order to fulfil kinematic admissibility (i.e., volume conservation) this drop in height results in soil crossing velocity discontinuity VD1, where its velocity changes from ‫ݒ‬ ଶ to ‫ݒ‬ ଵ and the soil is subjected to the plastic shear deformation ߛ ଵ = tan ିଵ ߯ ଵ = ߱ equal to the angle between the D r a f t 10 two velocities. At every displacement increment a new portion of the sheared material crosses VD1 and contributes to the evolving smeared shear zone. In contrast, no material crosses the velocity discontinuity VD2, where additional shearing takes place along a thin shear band. The plastic shear strain in VD2 can be determined using equation (1):

Work equation
The mechanism in Fig. 4 will, in general, accelerate due to softening of the material, with kinetic energy rate calculated as the difference between the rates of external work and dissipation: The weights of Blocks 1 and 2 have to be evaluated in order to calculate the change in their potential energy (work increments of gravity forces). The initial volumes of Blocks 1 and 2 are The volume of the sheared material which has crossed the velocity discontinuity VD1 from Block 2 to Block 1 is and is equal to the area of the grey smeared shear zone in Fig. 4.
For submerged slopes only the effective soil weight has to be taken into account for calculating the external work of the gravity. The rate of external work of the gravity acting on Blocks 1 and 2 and the reaction force ܳ ௨ acting on Block 1 is The dissipation rate along the slip surface is calculated from the velocity of Block 1 and the resistance ߬ on the slip surface.
The dissipation rate along velocity discontinuity between Block 2 and Block 3 (VD2 in Fig. 4) is calculated by integrating the product of the relative velocity ‫ݒ‬ ଶ and the corresponding shear resistance ‫ݏ‬ ௨ ൫ߛ ଶ ൯ along the length of the VD2, which changes depending on the drop in height of Block 2.
The dissipation rate in the shear zone between Block 1 and Block 2 (VD1 in Fig. 4) is obtained using equation (8): Substituting equations (9)-(12) into (6), the following expression for the kinetic energy rate is

Optimizing the mechanism geometry
Having formulated the terms of the work equation, the free parameter of the mechanism, the angle ߱, has to be optimized. The most probable geometry of the mechanism is found by following the minimum energy principle. This optimization is done for the initial failure by formulating the kinetic energy rate for zero displacement ‫ݑ‬ ଵ = 0. In this case, from equations (4), (5), (8) and (13) it follows that ‫ݑ‬ ௗ = 0, ܸ ௦௨ = 0, ߯ ଶ = 0 resulting in where ‫ݏ‬ ௨ = ‫ݏ‬ ௨ ሺ߯ ଶ = 0ሻ.
Following the minimum energy principle, the optimum angle ߱ ௧ is found from equation For perfectly plastic soil behaviour ‫ݏ‬ ௨ ሺ߯ሻ = ‫ݏ‬ ௨ = ‫.ݐݏ݊ܿ‬ and horizontal ground surface ߙ = 0, the above equation gives ߯ ௧ = 1 (i.e., the optimum angle ߱ ௧ = ߨ 4 ⁄ ), as in the classical undrained active failure mechanism. In the case of strain softening and mild slope inclinations ߙ, equation (15) gives ߯ ௧ ≅ 1, i.e., the shear strain along the velocity discontinuity between Block 1 and Block 2 (VD1 in Fig. 4) is approximately 100%. Such strains would most likely cause the full softening of the material: ‫ݏ‬ ௨ ሺ߯ ≅ 1ሻ = ‫ݏ‬ ௨ , so that the optimal geometry of the mechanism can well be approximated as Once the geometry of the initial failure mechanism has been defined, criteria for initial failure at the top of the sliding slab are formulated.

Conditions for the initial failure
Following the principles of limit analysis, the condition for initial (static) failure is determined by the energy balance equation: The initial failure will take place at the intersection of velocity discontinuities VD1 and VD2 in Fig. 4, where the slope parallel force ܳ ௨ in the sliding slab reaches its minimum and where we set the origin of coordinate ‫.ݔ‬ For optimized geometry and ‫ݔ‬ = 0, applying condition (17) to equation (14) gives the limiting slope parallel force in the sliding slab that would lead to the formation of the initial failure mechanism in the slab due to unloading: For perfectly plastic soil behaviour ‫ݏ‬ ௨ ሺ߯ሻ = ‫ݏ‬ ௨ = ‫ݐݏ݊ܿ‬ and horizontal ground surface ߙ = 0, equation (15) gives ߯ ௧ = 1, so that the limiting force (18) reduces to ܳ ௨ = ఊ ᇲ మ ଶ − ‫ݏ2‬ ௨ ℎ as in the classical undrained active failure mechanism.
Once the force ܳ ௨ is reached in the uphill part of the slab, the initial failure mechanism (Fig. 4) will form. Continuing downhill displacement of the slab will cause evolution of this mechanism which will result in the further unloading of the stable part of the slope.
D r a f t 14

Evolution of the initial failure mechanism at the top of the sliding slab
Once the force ܳ ௨ is reached in the uphill part of the slab, the initial failure mechanism (Fig. 4) will form. Continuing downhill displacement of the slab will cause further evolution of this mechanism. The direction and location of the velocity discontinuities VD1 and VD2 are not likely to deviate from the angle ߱ ௧ of the initial failure mechanism due to the material softening and the resulting shear strain localization. The displacements of the mechanism will result in the further unloading of the stable part of the slope and, potentially, propagation of the failure into this stable part. In the following section this unloading is quantified by calculating the slope parallel force ܳ , which is acting between Blocks 2 and 3 during the evolution of the initial unloading mechanism ( Fig. 5). The rate of kinetic energy for the evolving mechanism has been calculated from the energy balance equation (13) as a difference between the external work and internal dissipation rates. It can be also calculated by differentiating the total kinetic energy of the two Blocks 1 and 2 with respect to time: where the density of soil ߩ is introduced. Substituting expression (19) into the energy balance equation (13) provides the expression relating the slab force ܳ ௨ to the acceleration ‫ݒ‬ሶ ଵ and velocity Formulating equation of motion in the slope parallel direction includes inertia forces ‫ܫ‬ ଵ and ‫ܫ‬ ଶ , which can be found from time derivatives of the momentum of Blocks 1 and 2, respectively: Equation of motion in slope parallel direction is given by which after substitution of the slab force ܳ ௨ from equation (20) and inertia forces from equation (21) provides the expression for the force ܳ acting on the stable part of the slope as a function of the drop of the seabed surface ‫ݑ‬ ௗ , the velocity of the mechanism ‫ݒ‬ ଵ , slope geometry ߙ and shear strength of the soil ‫ݏ‬ ௨ : where This unloading force ܳ ሺ‫ݑ‬ ௗ , ‫ݒ‬ ଵ ሻ does not depend on the slab force ܳ ௨ , justifying the application of the free body cut of the slab in the above analysis. Knowing the evolution of this force is essential for formulating conditions at which the initially stable upper part of the slope (Block 3) becomes unstable, potentially leading to retrogressive failure.

Retrogressive failure mechanism
As the sliding slab continues to move downslope, the evolution of the initial failure mechanism in eventually causing the formation of the next failure mechanism, upslope from the initial one ( Fig.   6). This recurrent formation of failure mechanisms explains the process behind the retrogressive failures observed in spreads. The same large deformation KEB approach as for the initial failure mechanism can be applied to subsequent mechanisms. The mechanism in Fig. 6 is very similar to the initial failure mechanism in Fig. 4, and will fail when the limiting force ܳ ௨ from equation (18) is approached at the intersection of velocity discontinuities VD3 and VD4. This will happen when the force ܳ ሺ‫ݑ‬ ௗ , ‫ݒ‬ ଵ ሻ from equation (23), which supports the stable part of the slope, drops to the critical value ܳ , determined from the equilibrium of Block 3 (and a portion of Block 4 directly above it), provided the shear strength on the SPSS is known: Substituting the limiting force ܳ ௨ from equation (18) into equation (25) gives the critical value for the unloading force ܳ indicating initiation of the retrogressive failure above the sliding slab: where ߬ = ℎߛ ᇱ ‫݊݅ݏ‬ ߙ. The critical drop of the seabed level ‫ݑ‬ ௗ causing the retrogressive failure can then be found from equation (23) for ܳ ሺ‫ݑ‬ ௗ , ‫ݒ‬ ଵ ሻ = ܳ : where parameters a and b are given by equations (24), and Note that for perfectly plastic material ‫ݏ(‬ ௨ = ‫ݏ‬ ௨ = ‫,)ݐݏ݊ܿ‬ quasi-static unloading ‫ݒ(‬ ଵ = 0) and average shear strength on SPSS equal to the gravitational shear stress (݀ = 0), equation (27) gives ‫ݑ‬ ௗ = 0, i.e., no drop in the seabed surface is required for the continuing development of spreading failure.

Spreading failure criterion with progressive shear band propagation
From equation (26) it follows, that the critical load for initiating the retrogressive failure is dependent on the shear resistance on the SPSS. If Blocks 2, 3, and 4 are assumed to be rigid, D r a f t 18 limiting equilibrium considerations suggest that the failure would be initiated when the force ܳ ሺ‫ݑ‬ ௗ , ‫ݒ‬ ଵ ሻ from equation (23) drops to the value determined by the peak strength ߬ along the slip surface segment ‫ܮ‬ ଷ : The corresponding critical drop of the seabed level ‫ݑ‬ ௗ causing retrogressive failure can be found from equation (27) by substituting there ݀ = ൫߬ − ߬ ൯‫ܮ‬ ଷ .
In reality, however, the pre-failure unloading of the stable part of the slope results in its elastic deformation, which can lead to the progressive uphill propagation of a shear band along the SPSS with subsequent reduction of the shear resistance from its peak value ߬ to the residual value ߬ along the length of the propagated shear band ‫ܮ‬ ௌ௨ . This results in higher values of the critical force ܳ than predicted by the limiting equilibrium: i.e., for the initiation of retrogressive failure less unloading is required. Following Puzrin et al.
(2017) the propagation length of such a shear band in the stable zone can be quantified by means of the energy balance approach (Palmer and Rice, 1973): where ܳ is the initial slope parallel force in the slope, defined using the at rest earth pressure coefficient ‫ܭ‬ , ‫ܮ‬ is the critical length for catastrophic SBP, r is the gravitational shear stress D r a f t 19 ratio, E' is the plane strain unloading modulus and ߜ ̅ is the characteristic displacement (Fig. 3b) (23): where Parameters a and b are given by equations (24), ܳ ா by equation (29), r and ‫ܮ‬ by equations (31).
Note, that this approach assumes the process zone length to be small compared to the depth of the SPSS. In case the assumption of the small process zone length should not hold, an alternative process zone approach for assessing the shear resistance along the SPSS should be used (e.g., Puzrin and Germanovich 2005;Zhang et al. 2015).
D r a f t 20

Discussion of the spreading failure criterion
Consider a slope ߙ with the thickness of the sliding layer ℎ and the effective specific weight ߛ′.
The uniform peak undrained shear strength is assumed to be equal to that in the middle of the layer, calculated as a fraction ݇ of the corresponding vertical effective stress (residual strength is a fraction of the peak one, calculated using the sensitivity ‫:)ݏ‬ Localized failure along the SPSS requires that its shear strength is lower (a fraction Γ < 1) than the shear strength of the adjacent material in the sliding body: The softening behaviour of the material in the sliding body is assumed to be exponential (Einav and Randolph 2005), with the characteristic shear strain ߯ , The characteristic displacement of the SPSS and the stiffness of the sliding body are varied in a parametric study. The values of the parameters used in the parametric study presented below are shown in Table 1. These values lie within the range typical for Azeri submarine landslides investigated by Gray et al (2015) and Puzrin et al (2016). The thickness of the velocity discontinuity within the sliding body ݀ ௌ,ଶ has been chosen of the order of magnitude of 1 mm, typical for the fine grain sediments.
Insert Table 1 here.
For quasi-static conditions and the shear stress ratio ‫ݎ‬ = 0 ሺ߬ = ߛ ᇱ ℎ‫ߙ݊݅ݏ‬ሻ, the evolution of the unloading force ܳ from equation (23) and Fig. 6 as a function of the drop of the seafloor ‫ݑ‬ ௗ is presented in a normalized form as a thin solid black curve with dot markers in Fig. 7a. The initial D r a f t 21 failure forms at ‫ݑ‬ ௗ = 0, when the force in the slab at the intersection of VD1 and VD2 drops to ܳ ௨ defined from equation (18). The corresponding initial unloading force ܳ is found from equations (23) and (24) for ‫ݑ‬ ௗ = 0 and ‫ݏ‬ ௨ ሺ߯ ଶ ሻ = ‫ݏ‬ ௨ . As ‫ݑ‬ ௗ starts to increase, the unloading force ܳ first also increases due to the rapid softening ‫ݏ‬ ௨ ሺ߯ ଶ ሻ → ‫ݏ‬ ௨ along the velocity discontinuity VD2, decreasing coefficient b in equations (23) and (24)

Failure at the bottom of the sliding slab
In this section the KEB approach is used to formulate the conditions for initial failure at the bottom of the sliding slab, the subsequent evolution of the force in the lower part of the slope as a function of the state of plastic deformation and the criteria for the ploughing/runout failure.

Initial loading failure mechanism
The large deformation loading mechanism (Fig. 8) and the corresponding heave of Block 2 is With increasing deformation, the upper left corner of Block 2 will be overhanging and will fail at a certain stage of mechanism deformation. However, for assessing geometry and conditions for initial failure (i.e., small deformations), the potential failure of the corner can be neglected. For subsequent evolution of the mechanism at large deformations the failure of the corner will be included in the analysis.
The work equation can be formulated in the same way as for the unloading mechanism (6) and the initial volumes of Blocks 1 and 2 are given in equations (7). The volume of the material crossing the velocity discontinuity VD1 is ܸ ሶ ௦ = ℎ‫ݒ‬ ଵ ; ܸ ௦ = න ℎ‫ݒ‬ ଵ ‫ݐ݀‬ = ℎ‫ݑ‬ ଵ (38) and is equal to the area of the grey smeared shear zone in Fig. 8. The rate of external work with respect to the gravity acting on Blocks 1 and 2 and the driving force ܳ is given by The dissipation rates along the velocity discontinuities are formulated in analogy to the unloading mechanism ‫ܦ‬ ሶ ௌௌௌ + ‫ܦ‬ ሶ ଶ + ‫ܦ‬ ሶ ଵ = න ߬ሺ‫ݔ‬ሻ‫ݒ‬ ଵ ‫ݔ݀‬ The optimization of the mechanism geometry results in the optimal angle ߱ ௧ which is the same as for the unloading mechanism given by equation (15). Using condition (17) and in analogy to the unloading failure mechanism, the threshold for initial failure at the bottom of the sliding slab can be expressed by the limiting loading failure force

Evolution of the force acting on the stable part of the slope
In order to describe the evolution of the slope parallel force ܳ , which is acting on the stable part of the slope (i.e., Block 3), it is essential to include large deformations and geometric nonlinearities into the analysis. Similar to the unloading mechanism, the material crossing velocity discontinuity VD1 and increasing the volume of Block 2 has to be taken into account. The failure of the soil at the overhanging corner of Block 2 and the corresponding loss of volume of Block 2 should also be taken into account (Fig. 9). In general, the breaking of soil from the overhanging corner of Block 2 can be assessed by means of upper bound limit analysis, leading to a mechanism with an increasing number of failed blocks on the seafloor. This mechanism, however, can be drastically simplified by replacing the breaking of discrete blocks with continuous shearing of the soil, which crosses velocity discontinuity VDR and forms the new Block R. It can be shown that a good fit between the general and simplified mechanisms can be achieved by choosing a shearing angle as ߱ = ߱, resulting in the kinematic relation ‫ݒ‬ = ‫ݒ‬ ଵ , and the volume of the sheared material crossing velocity discontinuity VDR is given by As the slab keeps sliding and the new material is sheared into Block R, it moves over the top of the stable Block 3 with velocity ‫ݒ‬ = ‫ݒ‬ ଵ , overcoming resistance ‫ݏ‬ on the interface between blocks R and 3. This resistance, however, is likely to be very low due to the intercalation of the seawater between the original seabed and the overflowing sheared material. In the following we assume The rate of kinetic energy for the evolving mechanism can be calculated from the energy balance equation as a difference between the external work rates and internal dissipation rates The same rate of kinetic energy can also be calculated by differentiating the total kinetic energy of the three Blocks 1, 2 and R with respect to time: Substituting expressions (44)-(46) into the energy balance equation (43) provides the expression relating the slab force ܳ to the acceleration ‫ݒ‬ሶ ଵ and velocity ‫ݒ‬ ଵ of the mechanism: Formulating equation of motion in the slope parallel direction requires inertia forces ‫ܫ‬ ଵ , ‫ܫ‬ ଶ , and ‫ܫ‬ , which can be found from time derivatives of the momentum of Blocks 1, 2 and R, respectively: The equation of motion in slope parallel direction is then given by which after substitution of the slab force ܳ from equation (47) and inertia forces from equations (48) provides the expression for the force ܳ acting on the stable part of the slope as a function of the heave of the seabed surface ‫ݑ‬ , the velocity of the mechanism ‫ݒ‬ ଵ , slope geometry ߙ and shear strength of the soil ‫ݏ‬ ௨ :

Ploughing failure criterion
The mechanism in Fig 11 is very similar to the initial failure mechanism in Fig. 8, and will fail when the force acting on Block 3 will reach its critical value ܳ . This critical force is larger than the one for the first failure mechanism in Fig. 8, because by the time the second failure mechanism develops, there is an overflow of material from Block 2 over Block 4, increasing its volume: In addition, because of this overflow, the length of the velocity discontinuity VD3 is larger than that of VD1 in the initial mechanism, resulting in a higher rate of the volume crossing this discontinuity: D r a f t 28 Following the principles of limit analysis, the condition for initial (static) failure of the second mechanism is determined by the energy balance equation for Blocks 3 and 4: The rate of external work with respect to the gravity acting on Blocks 3 and 4 and the critical driving force ܳ is given by The dissipation rates along the velocity discontinuities are formulated in analogy to the first failure mechanism: Substituting equations (54) and (55) into (53) gives the value of the critical driving force ܳ causing the initiation of the second (and all subsequent) ploughing failure mechanisms: where ߬ = ℎߛ ᇱ ‫݊݅ݏ‬ ߙ. This second mechanism will develop when the force ܳ ሺ‫ݑ‬ , ‫ݒ‬ ଵ ሻ from equation (50), which loads the stable part of the slope, reaches the critical value ܳ from equation (56). This takes place, when the heave of the seabed level ‫ݑ‬ reaches its critical value: Note, that similar to spreading, in the case of a perfectly plastic material ‫ݏ(‬ ௨ = ‫ݏ‬ ௨ = ‫,)ݐݏ݊ܿ‬ quasi-static loading ‫ݒ(‬ ଵ = 0) and average shear strength on the SPSS equal to the gravitational shear stress (݀ = 0), equation (57) gives ‫ݑ‬ = 0, i.e., no heave in the seabed surface is required for the continuing ploughing. From equation (56) it can be seen, that the failure load causing ploughing depends on the shear resistance on the SPSS. If the Blocks 2, 3, and 4 are assumed to be rigid plastic, limiting equilibrium considerations suggest that the failure would be initiated when the force ܳ ሺ‫ݑ‬ , ‫ݒ‬ ଵ ሻ from equation (50) drops to the value determined by the peak strength ߬ along the slip surface segment ‫ܮ‬ ଷ : The critical heave of the seabed level ‫ݑ‬ causing ploughing is then defined from equation (57), In reality, however, the pre-failure loading of the stable part of the slope results in its elastic deformation, which can lead to progressive downhill propagation of a shear band along the SPSS D r a f t 30 with subsequent reduction of the shear resistance from its peak value ߬ to the residual value ߬ along the length of the propagated shear band ‫ܮ‬ ௌ . This results in lower values of the critical force ܳ than predicted by the limiting equilibrium: i.e., for the initiation of the ploughing smaller loads are required. Following (Puzrin, 2016) the propagation length of such a shear band in the stable zone can be quantified by means of the energy balance approach (Palmer and Rice, 1973): where ܳ is the initial slope parallel force in the slope. The critical length for catastrophic SBP ‫ܮ‬ and the gravitational shear stress ratio r are given by equations (31). In case ‫ܮ‬ ௌ > ‫ܮ‬ ଷ = ℎ ఞ మ ାଵ ఞ , the whole base ‫ܮ‬ ଷ of the Block 3 fully softens to ߬ and the corresponding critical heave of the seabed level ‫ݑ‬ causing ploughing can be found from equation (57) by substituting there ݀ = ൫߬ − ߬ ൯‫ܮ‬ ଷ . In contrast, for ‫ܮ‬ ௌ < ‫ܮ‬ ଷ , only the portion ‫ܮ‬ ௌ of the base ‫ܮ‬ ଷ fully softens, and a larger heave of the seabed ‫ݑ‬ is needed in order to overcome the peak strength ߬ of the remaining portion ‫ܮ‬ ଷ − ‫ܮ‬ ௌ of the base. This heave can be found from the condition ܳ ሺ‫ݑ‬ , where the critical force ܳ ா also depends on the loading force ܳ ሺ‫ݑ‬ , ‫ݒ‬ ଵ ሻ. After resolving the ploughing condition ܳ ሺ‫ݑ‬ ௗ , ‫ݒ‬ ଵ ሻ = ܳ ா with respect to ܳ ሺ‫ݑ‬ , ‫ݒ‬ ଵ ሻ, the critical heave of the seabed surface can be found from the following quadratic equation: D r a f t 31 derived using equations (50) and (59) -(61). Parameters r and ‫ܮ‬ are given by equations (31).
Note, that this approach assumes the process zone length to be small compared to the depth of the SPSS.

Discussion on ploughing criterion
Similar to the spreading case in Fig. 7, a limited parametric study (with the same set of parameters from Table 1) is used here to discuss the ploughing criteria. Compared to the spreading case in Fig.   7, the range of the combination of SBP parameters ‫ߜܧ‬ ̅ /ሺߛ ᇱ ℎ ଶ ሻ in Fig. 12 has been shifted towards higher values, in order to obtain the complete range of the normalized critical heaves. For quasistatic conditions and the shear stress ratio ‫ݎ‬ = 0 ሺ߬ = ߛ ᇱ ℎ‫ߙ݊݅ݏ‬ሻ, the evolution of the loading force ܳ from equation (50) and Fig. 10 is presented in a normalized form as a function of the heave of the seafloor ‫ݑ‬ (thin solid black line with dot markers in Fig. 12a). The initial failure forms at ‫ݑ‬ = 0, when the force in the slab at the intersection of VD1 and VD2 reaches ܳ defined by equation (41). As ‫ݑ‬ starts to increase, the loading force ܳ first decreases due to the rapid softening ‫ݏ‬ ௨ ሺ߯ ଶ ሻ → ‫ݏ‬ ௨ along the velocity discontinuity VD2. Further deformation leads to increase of the normalized loading force caused by the geometric hardening of the sliding body according to equation (50). Similar to spreading, for a dynamic case with a constant velocity of Block 1 of ‫ݒ‬ ଵ = ‫,ݏ/݉2‬ inertia of the failing blocks lead to an increase of the loading force (dashed line), but this dynamic effect is sufficiently small to justify a quasi-static approximation of the ploughing failure criterion.
The left and right thin dashed curves in Fig. 12a represent the two extreme cases of the critical force ܳ ா from equation (60), which are required to initiate ploughing in the sliding body. The curved shape results here from the fact that in contrast to spreading, not only the loading force ܳ but also its critical failure limit ܳ ா are functions of the heave ‫ݑ‬ , which is a result of the material being deposited on top of the subsequent failure mechanism. The right curve corresponds to the D r a f t 32 case of ‫ܮ‬ ௌ = 0, i.e., the peak strength ߬ is mobilized along the entire base ‫ܮ‬ ଷ of Block 3 in Fig   11. The left curve presents the fully softened case of ‫ܮ‬ ௌ = ‫ܮ‬ ଷ , with only residual strength ߬ mobilized along ‫ܮ‬ ଷ . Intersection of the thin solid black line (the loading force ܳ ) with the left curve (point A) represents the value of the minimal critical heave ‫ݑ‬ for the case of the shear band propagating far into the stable part of the slope. The intersection of the thin solid black line with the right curve (point B) represents the value of the minimal critical heave ‫ݑ‬ , that will cause ploughing even when no SBP takes place. In between these two extreme cases, the shear band propagates over a portion of the base of Block 3 ሺ‫ܮ‬ ௌ < ‫ܮ‬ ଷ ሻ and the critical heave ‫ݑ‬ is found from the ploughing criterion in equation (62) as a function of the SBP parameters ‫ߜܧ‬ ̅ /ሺߛ ᇱ ℎ ଶ ሻ. In While in a curved slope it is likely that the initial failure of the sliding body occurs at ‫ݎ‬ = 0, further ploughing may progress into portions of the slope with lower or higher inclination. The corresponding required critical heave for ploughing is shown in Fig. 12b as a function of the shear stress ratio ‫.ݎ‬ For each value of r, the lower solid black curve presents the fully softened case ‫ܮ(‬ ௌ = ‫ܮ‬ ଷ ), corresponding to the point A in Fig. 12a Example of a post-failure evolution in a curved slope

Model parameters
The evolution of a submarine landslide using a numerical incremental integration of the developed KEB approach is illustrated by the following example. Consider the skew-symmetric S-shaped slope in Fig. 13 with geometry built from a linear part (inclination ߙ and length ‫ܮ2‬ ) in the middle with a smooth transition to milder slopes up-and downhill, described by the following exponential functions: The geometric parameters are given in Table 2.The shear strengths on the SPSS and in the sliding body have been adopted from  Fig. 13.
Insert Table 2 here.

Incremental integration of the dynamic landslide evolution
As the first step, we introduce an initial shear band on the SPSS in the unstable region, where the material has fully softened: Stage 1 in Fig. 13. This leads to elastic unloading at the top and loading at the bottom end of the shear band. The force in the sliding body is compared to the criteria for D r a f t 34 initial unloading and loading failure in equations (18) and (41), respectively. In the presented example the unloading failure criterion (18)  Stage 2 in Fig. 13. However, even before this passive failure at the bottom takes place, the elastic deformation of the slab leads to further unloading and a small drop in height in the active failure mechanism at the top, related by equation (23). Now that the global failure up-and downhill takes place, the continuing downslope movement of the slab causes an additional drop of the seabed surface in the active mechanism and a heave in the passive mechanism, leading, respectively, to further unloading and loading described by equations (23) and (50). After each displacement increment the conditions (32) for spreading at the top and (60) for ploughing at the bottom are checked, and if fulfilled, the next failure mechanism is formed. In the presented example, first, in parallel with the uphill retrogressive spreading, the uphill (retrogressive) ploughing occurs: Stage 3 in Fig. 13. (Note, that the criterion for retrogressive ploughing is slightly different to that in equation (60). Later, also downhill ploughing takes place and, due to the significant downhill slab displacements, the retrogressive spreading blocks in the upper part of the slope start moving downslope: Stage 4 in Fig. 13. The total kinetic energy of the system is calculated by summing up incrementally the kinetic energy of the moving slab together with the kinetic energy of all the spreading and ploughing mechanisms, calculated using equations (13) and (43)

Summary and conclusions
The paper applies the kinematic energy balance (KEB) approach based on the kinematic method of plasticity theory to the large deformation problem of initiation and propagation of the spreading and ploughing failure, at the top and bottom edges of a failed submarine slope, respectively. The models account for the phenomenon of a progressive propagation of a slope parallel slip surface, which is also quantified using the energy balance approach. The investigation of the evolution of the forces and post-failure geometry in the sliding layer allows for formulating spreading and ploughing criteria both in terms of the critical lateral forces and the critical drop (for spreading) and heave (for ploughing) of the seabed surface. These criteria appear to be strongly dependent not only on the shear strength on the slip surface and in the sliding body, but also on the stiffness of the sliding body and the strain softening in the process zone of the propagating slip surface.
In contrast to existing approximate ad-hoc analytical and numerical solutions, the proposed approach provides rigorous theoretical basis for spreading and ploughing criteria. Furthermore, it has been demonstrated that the incremental energy balance equations, which were derived to define the above criteria, can also be used for deriving a comprehensive dynamic solution of the postfailure landslide evolution. Incremental integration of the kinetic energy in time allows for monitoring the development of normal forces in the sliding slab leading to the initiation of new kinematic failure mechanisms. Once a new mechanism forms, its deformation can be modelled using the same energy balance equations, resulting in the final post-failure geomorphology of the failed slope. Note, however, that because the mechanisms used in the large-displacement kinematic analysis, although admissible, are not necessarily true, the displacements calculated from such kinematics may become progressively inaccurate, even if the boundaries are continuously updated D r a f t 36 during the process. Therefore, the proposed mechanisms, while providing a more rigorous plasticity approach to the problem, require further validation against the field observations.
Although the overflow of the sheared material over the stable slope below the failed slab has been incorporated into the ploughing mechanism, the limited scope of this paper did not allow for presenting derivations for the corresponding run-out criteria. These criteria will be presented in a subsequent publication that will complete the formulation of the KEB framework for the postfailure landslide evolution and enable its application to case studies.     D r a f t D r a f t Table 2. Geometric parameters for the example of submarine landslide evolution.