Frictional heating lubrication for submarine landslide

Given a large block of fine-grained sediment that overlies a 1° continental slope, if the entire sediment deposit is shaken by an earthquake to slide down the slope with an initial velocity (referred to as “submarine landslide” hereinafter), our block model showed that the whole deposit can continuously slide downslope with the “help” of basal frictional heating. In theory, basal Coulomb friction can generate a thermal internal energy in the basal shear zone of landslides (called “frictional heating”), the increasing temperature activating two mechanisms. One mainly decreases the Terzaghi effective normal stress of the landslides and the other decreases the Coulomb friction coefficient. Combining these two mechanisms can effectively decrease the basal frictional resistance of the landslides increasing the mobility of the landslides on gentle slopes (entitled as “frictional heating lubrication”). As shown in our calculations, frictional heating increases with the increase of the initial downslope velocity but decreases with the increase of the shear zone thickness. Article history: Received 15 April 2016 Revised 21 March 2017 Accepted 22 March 2017


INTRODUCTION
Sediments derived from land and from the continental shelf are often deposited on the upper continental slope.The former is mainly transported by rivers whereas the latter by ocean currents or storms.Once released, the sediments usually flow in a huge volume down the slope (Masson et al. 2006).In this study, we define such kind of massive sediment movement as "submarine landslide", which often can travel a long distance even on gentle slopes.To date, it is widely accepted that the occurrence of most submarine landslides is not mainly influenced by slope gradients.In fact, many studies have shown that submarine landslides often occur on slopes the gradients of which are as low as 1° (Hampton 1972;Embley 1976;Hampton et al. 1996;Vorren et al. 1998;McAdoo et al. 2000;Huhnerbach et al. 2004).Submarine landslides are one of the effective processes that can transfer tremendous sediments across the continental slope to the deep ocean.The largest submarine landslides can involve thousands of km 3 of material, two to three orders of magnitude larger than any terrestrial landslides (see Table 1, Hampton et al. 1996).For example, the Storegga slide in-volved about 3000 km 3 of sediment, affected 95000 km 2 of the Norwegian slope and basin, and had a run-out distance of around 800 km (Haflidason et al. 2004).This large area is about 20% bigger than Scotland and the distance is close to the length of mainland Britain.Therefore, because of large volumes and long traveling distances, submarine landslides not only threaten to destroy coastal settlements and offshore oil installations but can also destruct submarine telephone cables far from the shore.More severely, submarine landslides can generate tsunamis to cause tremendous casualties.For instance, the 1998 tsunami striking Papua New Guinea killed 2200 people there (Satake and Tanioka 2003).
Many studies have suggested possible factors that can trigger submarine landslides, e.g., earthquakes, hurricanes or cyclic loading, oversteepening of slopes, overpressure, gas hydrate dissociation, sea-level change, and volcanic activity (Prior and Coleman 1982;Weaver and Kuijpers 1983;Moore et al. 1989;Assier-Rzadkiewicz et al. 2000;Sultan et al. 2004a, b;Fine et al. 2005).However, it is still a question why submarine landslides can occur on slopes as low as 1° and travel a long distance over such mild slopes.Recently, a model called "hydroplaning lubrication theory (HLT)" has been proposed to explain why submarine landslides Terr. Atmos. Ocean. Sci., Vol. 29, No. 1, 87-103, February 2018 comprised of finely comminuted clay can travel far on gentle slopes (Mohrig et al. 1998;De Blasio et al. 2004), the HLT model briefly explained as follows.When submarine landslides are in motion, according to the non-Newtonian fluid mechanics, the vertical depth of the landslides can be divided into two layers: the upper one is an "unsheared plug zone" and the lower one a "shear zone".The HLT, neglecting the shear zone, assumes that the entire depth is a plug zone (De Blasio et al. 2004).A more peculiar assumption is that when a landslide reaches a critical velocity, water at the front stagnation point with the highest pressure will intrude the bottom of the landslide and penetrate underneath until it reaches the back of the landslide.Based on this assumption, there will be a water layer existing between the bottom face of the landslide and the seabed surface and acting as a lubricant to reduce the frictional resistance of the landslide.However, this assumption could be a weakness of the theory.For example, given a large landslide with a depth 200 m moving on a gentle slope, according to Eq. ( 3) (De Blasio et al. 2004), the critical velocity for this landslide to activate the hydroplaning effect is about 18 m s -1 .When water intrudes the bottom of the landslide at this high velocity, it is more likely that the water will scour the seabed and entrain the basal sediments to form a solids-liquid (twophase) shear flow there, instead of a water layer.Moreover, the HLT model does not explain how submarine landslides can overcome the basal frictional resistance and speed up from a "static state" to reaching a high critical velocity (e.g., 18 m s -1 ).Up till today, no data of large-scale landslides have been found to support the HLT yet.In fact, hydroplaning of natural debris flows has never been directly observed.
Our motivation in this article is to seek a model that can explain why submarine landslides can move on gentle slopes.Our model is different from the HLT.We do not assume a water lubricating layer.In our model, the entire landslide depth is comprised of two parts: an upper "insignificant shear zone" and a basal "shear zone".In this basal shear zone, water and sediments get mixed as a solids-liquid (two-phase) flow, instead of a water lubricating layer (onephase flow).In this study, we show that the "frictional heating" in the basal shear zone is a key factor to the continuous movements of submarine landslides on gentle slopes.We take account of three mechanisms (thermopressurization weakening, shear dilatancy, and thermal softening) in our model and study their impacts on the occurrence of submarine landslides on gentle slopes.
Thermopressurization weakening (TPW) mechanism has been studied for years in fault mechanics, particularly the gouge friction in the shear zone of faults (Segall and Rice 1995;Garagash and Rudnicki 2003;Rempel and Rice 2006;Rice 2006;Urata et al. 2013;Rice et al. 2014;Suzuki and Yamashita 2014).Such a highly localized shear layer has a thickness often on the order of centimeters, millimeters, or smaller, nested within the fault core (Rice 2006;Rice et al. 2014).It has been widely accepted that during the rapid shear of a fluid-saturated fault, the friction in the shear zone can cause an increase of temperature (called as "frictional heating"), which can induce a rise of excess pore pressure there.As a result, the Terzaghi effective normal stress (Skempton 1960) and thus the Coulomb frictional stress will decrease accordingly, leading to a fault slip.This process is referred to as thermopressurization weakening (TPW) mechanism, which has also been operated in the research of large terrestrial landslides (Voight and Faust 1982;Davis et al. 1990;Vardoulakis 2000Vardoulakis , 2002;;Goren andAharonov 2007, 2009;Veveakis et al. 2007Veveakis et al. , 2010;;Alevizos et al. 2010).In theory, the TPW mechanism acts as an "accelerator" that can speed up the movements of landslides, whereas the shear dilatancy mostly can behave as a "brake" to decelerate landslide movements (Iverson et al. 2000;Iverson 2005).In addition, we take account of another "accelerator" called as "thermal softening mechanism".When temperature is increased, for some kinds of clays, their coefficient of friction will decrease (Hicher 1974;Fig. 2 in Veveakis et al. 2010).This phenomenon is referred to as "thermal softening" (Vardoulakis 2000(Vardoulakis , 2002;;Veveakis et al. 2007Veveakis et al. , 2010;;Alevizos et al. 2010).To summarize, our model shows that the TPW and thermal softening mechanisms can effectively decrease the basal frictional resistance of submarine landslides.The former mainly increases the basal excess pore pressure to reduce the Terzaphi effective normal stress of landslides, whereas the latter decreases the coefficient of friction.These two mechanisms are activated by the basal frictional heating.
This article is organized as follows.In section 2, we introduce the governing equations.Our computational results are discussed in section 3, major conclusions presented in section 4.

MATHEMATICAL MODEL
Submarine landslides are often modeled as (1) an incompressible Newtonian viscous laminar flow (Jiang andLeblond 1992, 1994;Heinrich et al. 2000), (2) a non-Newtonian mud flow (Jiang and Leblond 1993;Assier-Rzadkiewicz et al. 1997), or (3) a sliding or slumping rigid block (Harbitz 1992;Grilli andWatts 1999, 2005;Hürlimann et al. 2000;Ward 2001;Grilli et al. 2002;Lynett andLiu 2002, 2005;Liu et al. 2005;Enet and Grilli 2007;Hornbach et al. 2007;De Blasio 2011;Wang et al. 2011;Gargani et al. 2014).As was described in De Blasio (2011, p. 319), it is often impossible to assess the rheology of submarine landslides.Not only are their material properties often inaccessible but their compositions also imprecisely known.Therefore, the first basic problem to simulate submarine landsides as a fluid is that the "stress-strain rate" constitutive relationships are unknown or unclear.Compared to the Navier-Stokes equations, block models are relatively simple and have been used widely, where the velocity of landslides is described by Newton's second law and only a few parameters are required.In the present study, our landslide model is basically the same as the "thermoporoelastic" slide block model which Goren andAharonov (2007, 2009) proposed.In addition, we incorporate the shear dilatancy and thermal softening mechanisms into our model to study the impacts of these mechanisms on the initiation of submarine landslides.We apply the concept of Goren and Aharonov to simulate submarine landslides as an intact slide block with most shear strains localized to a thin "zone" at the base of the block.This "shear strain localization zone", simply called as "shear zone" hereinafter, is not entirely an assumption.As explained below, this "shear zone" has been observed at the base of many subaerial landslides.This basal shear zone is a main part of our landslide model; as explained in section 2.3, the pore pressure generated by the basal frictional heating in the shear zone can significantly decrease the basal frictional resistance of landslides (called "TPW mechanism").As shown in our study here, the occurrence of submarine landslides on nearly flat seafloor is deeply affected by the existence of this basal shear zone.Understanding of this shear zone is helpful to compute the basal frictional resistance and acquire an insight into the triggering mechanism of submarine landslides.First of all, in section 2.1 we introduce the basal shear zone and present a literature review on "contracting behaviors" of strain localization zones, and subsequently, in sections 2.2 -2.4,we derive our landslide model including the governing equations for the basal shear zone.

Shear Zone Introduction
Shear zone (or called shear band) is an accepted concept in soil mechanics and particularly in slope analysis (Chowdhury et al. 2010, Chapter 8).It implies that within a sliding soil mass, most deformation is concentrated to a narrow "basal shear zone" as described earlier which separates the rest of the sliding mass where deformation may be insignificant.As shown in the observations of many subaerial landslides, the thickness of this "basal shear zone" is often very small ranging between mm-scale and decimeter-scale.For example, (1) an excavation of a failed clay slope in England artificially triggered by a fluid injection discovered that the basal shear zone was nested within a very thin disturbed zone (Cooper et al. 1998) ; (2) an investigation of a landslide in Hong Kong naturally triggered by a heavy rainfall found out that the basal shear zone of this landslide might be a few millimeters thick appearing at 2 -5 m depth below the slope surface (Wen andAydin 2003, 2004); (3) a four-year monitoring of a landslide in Italy, during which inclinometers were inserted into the landslide to measure the slide-soil displacement vs. the slope-normal depth, released that the landslide moved as a stiff body over a narrow zone (≈ decimeters thick) with most shear deformation localized to this thin zone (Picarelli et al. 2005); (4) a landslide experiment carried out on a natural slope, in which soil-strain probes were inserted into the sliding soil to measure the shear deformation, showed that the basal shear zone was only about 10 cm thick (Ochiai et al. 2004); (5) a field investigation of a sensitive clay landslide with a volume of 2.5 × 10 6 m 3 which occurred near British Columbia, Canada found out that the thickness of shear zone in this landslide could be narrowed down to a very thin "slip surface" (cf.Figs. 9, 10, or 11 in Geertsema et al. 2006); (6) two investigations of clay-slide deposits in Norway exhibited that the shear zone thickness had a range between millimeters and centimeters depending on the sediment properties and failure mechanisms (cf.Figs. 4 and 5 in Hansen et al. 2007;Fig. 7e in Solberg et al. 2008); and (7) by installing inclinometers 50 m below the ground surface to measure the soil displacements at different slope-normal depths, a seven -year monitoring of a creeping, deep seated clayey landslide (38 m deep) in Italy displayed that even though the landslide movement is very slow, a shear zone with a thickness < 50 cm was clearly detected at the base of the landslide (Di Maio et al. 2013).
Next, we briefly describe the difference between "weak layer" and "shear zone".As revealed in some field investigations, downslope motion of submarine landslides might occur as a slip over a "weak layer" where the soil strength is very low or exceeded by gravitational driving force (e.g., Edwards et al. 1995;Gee et al. 2005;Blum et al. 2010).The existence of a "weak layer" may be a necessary condition for the initiation of submarine landslides.For submarine landslides, the first mention of "weak layer" might be found in Lewis (1971), but for subaerial landslides, it was invoked earlier by Kardos et al. (1944).Even for snow avalanches, "weak layer" was also well recognized (Heierli et al. 2008).
As far as we know, there are several types of "weak layers" observed in field investigations, for example, (1) sandy layers where gas hydrates dissociate (Sultan et al. 2004b); (2) consolidated clayey areas of high sedimentation rates (Dugan and Flemings 2000;Sultan et al. 2004a); (3) sandy silt layers that can liquefy during earthquakes (Kokusho and Kojima 2002); (4) rapid accumulation events, e.g., turbidity currents or debris flows, over a clay-rich seabed (Hansen et al. 2011); and (5) sensitive clay layers; i.e., clayey sediments with a "strain-softening" behaviour (Kvalstad et al. 2005;Dey et al. 2016).The "strain-softening" implies that just a small distortion can cause a tremendous deformation of soil (like a "collapsing" or "contracting" of soil texture) leading to a significant reduction of the porosity, which for undrained conditions will induce an increase of pore pressure there and therefore decrease the soil strength (Thakur 2011).Normally, just a small shear distortion can significantly reduce the soil strength of an extra sensitive clay.Due to this "strain-softening" property, the existence of sensitive clay layers is very possibly a main cause of many large-scale submarine landslides that occurred on nearly flat continental slopes, e.g., the Storegga slide (Bryn et al. 2005;Gauer et al. 2005;Kvalstad et al. 2005;Dan et al. 2007;Locat et al. 2014;Dey et al. 2016).For this reason, many scientists had proceeded to study sensitive clay.One of the main studies is to observe the "strain-localization-zone contracting" phenomenon in sensitive clays that is caused by the "strainsoftening" behaviour of the clays.As shown in the numerical study by Thakur (2011) and laboratory observation by Gylland et al. (2013), at the onset of shear strain localization, the localized zone in sensitive clay was a few millimeters thick, but this strain localization zone did not stop contracting until it was reduced to the minimum sizes.This "strain-localization-zone contracting" phenomenon not only proves that the clay has the "strain-softening" property but importantly, implies that the shear zone in sensitive-clay landslides can be very thin just like the example we cite above (Geertsema et al. 2006, see Figs. 9, 10, 11 therein).Recently, carrying out an experiment of granular shear flow, Van Der Elst et al. (2012) also discovered a "strain-localization-zone contracting" phenomenon that is caused by a granular rheology, not by the "strain-softening" behaviour.In earlier research (Mülhaus and Vardoulakis 1987), the strain localization zone observed for sand had a thickness on the order of tens times the mean particle size, the thickness probably ranging between millimeters and centimeters.But for kaolin clay, Vardoulakis (2000Vardoulakis ( , 2002) ) suggested the minimum localized thickness between microns and millimeters.To gain more information on shear behaviors of natural landslide materials or laboratory man-made soil, interested readers can be referred to the following microstructural analyses and references therein: Morgenstern and Tchalenko (1967); Foster and De (1971); Oda (1972); Lupini et al. (1981); Dewhurst et al. (1996); Anson and Hawkins (1999); Frost and Jang (2000); Wen andAydin (2003, 2004); and Haines et al. (2013).
The calculated results for the thickness of shear zone often vary with models (e.g., Platt et al. 2014;Chen and Rempel 2015).To date, it is still an open issue to predict the thickness of the basal shear zone in landslides because there are so many factors involved to affect the evolution of this shear zone.Therefore, in simulating the movements of large-scale subaerial landslides, researchers simply assume that the thickness of the basal shear zone is constant.For instance, in Veveakis et al. (2007) and Goren and Aharonov (2009) studying the 1963 Vaiont slide, the thickness was assumed as a constant equal to 0.1 and 0.161 m, respectively; in Goren and Aharonov (2007) simulating a long traveling landslide, the thickness was also assumed to be a constant (= 0.02 m); in Vardoulakis (2002), it was suggested that the basal shear zone had a fixed thickness between 0.0015 and 0.1 m.Accordingly, in our present research, we follow these studies to assume that (1) when a submarine landslide triggered to move downslope, shear deformation is instantly localized to a thin "zone" at the base of the landslide and (2) this basal shear zone has a constant thickness.To put more rigorously, in this preliminary research, our objective is to study "the initiation of submarine landslides affected by frictional heating in the case of constant basal shear zone thickness".For "non-constant" cases, interested readers can revisit our study by another paper.

Governing Equation for the Landslide Velocity
First of all, we assume that a parallelepiped block of marine sediment statically overlies a gentle slope (e.g., 1°).The dimensions of the block can reach several kilometers in length (L) and width (W) and approximately hundreds of meters in depth (D).Assume that at time (t) = 0, this sediment block was seismically shaken to slide downslope as an "intact" body.Based on Newton's second law, at time (t) > 0, the downslope velocity of this parallelepiped sediment block, denoted by u, is governed by the following equation, ( , ) where add c is the added mass coefficient, t D is is the density of this block that is a mixture of water and fine soil grains (called sediment), f t and s t represent the intrinsic density of the water and of the sediment, respectively, n is the porosity, z is the volume fraction of the fine soil grains (z = 1n), g = 9.8 m s -2 is the gravity, and i is the slope angle.sin g t i D and cos g t i D represent the slope-parallel and slope-normal components of the weight of this parallelepiped sediment block, respectively.Note that in these two components the hydrostatic buoyancy has been subtracted from the gross weight of the sediment block.The basal frictional resistance of this sliding sediment block is taken to be the product of the friction coefficient ( n) and the Terzaghi effective normal stress of the block, cos gD t i D p(z = 0, t), where p is the pore water pressure in excess of the hydrostatic pressure.The last term on the right-hand-side (RHS) of Eq. ( 1) accounts for the total water resistance caused by the form drag and skin friction, C D and C K being the drag and the skin-friction coefficients, respectively.
Basically, our block model has the following limits.(1) The thickness of the basal shear zone (h) has to be small because the strain rate (c o ) in the basal shear zone is approximated as linear (u/h).( 2) Water scouring effect has to be small and negligible because we assume that submarine landslides move on a gentle slope as an "intact" slide block.
(3) The slope surface has to be rigid because in our block model we do not consider sediment entrainments occurring between the bottom face of the slide block and the slope surface.Although having these limits, our landslide model can be applied to study a landslide-tsunami problem as described below.First, we all know that submarine landslides can generate large tsunamis to cause tremendous casualties.The damage degree is related to the height and velocity of the tsunamis, these two variables determined by the characteristics of the landslides, e.g., the velocity, dimensions (shape and volume), and subaqueous depth of the landslides.This is the traditional conclusion of landslide-tsunami studies.But if our present landslide model is coupled with a hydrodynamic model to compute the free surface waves generated by submarine landslides (e.g., one-way coupling), the results could be interesting.For instance, one can observe how the height and velocity of the surface waves are influenced by the basal -shear -zone thickness, frictional heating, thermopressurization effect, soil dilatancy, dissipation of excess pore pressure, etc., which has not been discussed yet.Particularly, one can see how the frictional heating lubricates submarine landslides and affects the velocity of the landslide-generated waves.

Governing Equation for the Excess Pore Pressure in the Basal Shear Zone
In this study, we consider that the excess pore pressure (p) in the basal shear zone is generated via several mechanisms such as thermopressurization, dilatancy, sedimentation, etc., the details discussed in this subsection.First of all, by assuming that water and fine sediments get mixed in the basal shear zone, the mass conservation of water in the pore space is expressed as (Rice et al. 2014) Where f t represents the intrinsic density of the water, n is the pore volume fraction (i.e., porosity), q f is the Darcy's flux of the water along the slope-normal coordinate (z), and h is the thickness of the basal shear zone.We follow Goren andAharonov (2007, 2009) to assume that all variables in the basal shear zone are functions of only time (t) and the coordinate z.Based on Darcy's law, q f is written as where l is the permeability, f h is the viscosity of the water, and p is the excess pore water pressure.Substituting Eqs.(3) into (2), one can rewrite Eq. (2) as According to the thermodynamics, the rate of change of the density f t could be expressed as the sum of the contributions from the changes of excess pore water pressure and temperature (T), where f b and f a are the compressibility and thermal expansivity of water, respectively.The first (or second) term on the RHS of Eq. ( 5) expresses that an increase of pore water pressure (or temperature) can increase (or decrease) the density f t .Substituting Eqs. ( 5) into (4), one can obtain the governing equation for the excess pore pressure p in the basal shear zone, The rate of change of temperature in Eq. ( 6) can be calculated using the following heat transfer equation, where c t denotes the sum of c n f f t and c s s t z, c f and c s are the specific-heat constants of the water and the sediment, respectively, x is the shear stress of the water-sediment mixture in the basal shear zone, and c o is the shear strain rate of the mixture.In Eq. ( 7), we neglect the heat convection and conduction because they are small within short landslide duration (Goren and Aharonov 2007).Equation ( 7) indicates that the temperature (T) in the basal shear zone is increased by the friction work xc o .This "temperature increasing" phenomenon is referred to as "frictional heating".In Eq. ( 7), the shear stress x has the same definition as was mentioned in Eq. ( 1), taken to be the product of the friction coefficient n and the Terzaghi effective normal stress cos gD p t i D -.The strain rate c o is assumed as u/h if the shear zone thickness (h) is small enough.
Here, we explain how we express the rate of change of porosity in Eq. ( 6) and the friction coefficient n.First of all, let us review an article by Boyer et al. (2011).In this article, the authors conducted ring-shear experiments of dense, non-Brownian hard spherical particles suspended in a viscous fluid which is sheared at a constant confining pressure.Boyer et al. (2011) did the experiments with different particle size, material, fluid viscosity, and initial number of particles.They obtained that all the suspensions can be well described by two constitutive relations: where x, eff v , n, and z are the shear stress, effective normal stress, friction coefficient, and solids volume fraction of the suspensions, and I v is the so-called "viscous number" defined as All the data obtained by Boyer et al. (2011) for different particle size, material, fluid viscosity, and number of particles show that n and z can be modeled as functions of I v , .
where n 1 = 0.32, n 2 = 0.70, and I 0 = 0.005 were chosen by Boyer et al. (2011) to agree with the experimental results of submarine granular flows (Cassar et al. 2005), and m z = 0.585 is the maximum value of z reached when I v approaches zero (quasi-static state).
Equations ( 8) and ( 9) are similar to the constitutive relations obtained in the case of dry dense granular flows (GDR MiDi 2004;Da Cruz et al. 2005;Jop et al. 2006;Forterre and Pouliquen 2008), which had been proved very successful in describing the rheology of dry dense granular flows in a wide range of flow configurations (Forterre and Pouliquen 2008).Boyer et al. (2011) demonstrated that the rheology of non-Brownian dense spherical particles suspended in viscous fluids is like the rheology of dry dense granular flows, both rheologies completely determined by a friction law and a volume-fraction law.The friction law which Boyer et al. (2011) proposed shown in our Eq.( 8) exhibits that the shear stress ( x) of the dense suspensions increases with the increase of the effective normal stress ( eff v ) or the friction coefficient ( n) that is expressed in Eq. ( 11) as an increasing function of the "viscous number" (I v ).On the RHS of Eq. ( 11), the sum of the first two terms is a contribution from the rubbing contacts among dense granular solids, whereas the last two terms represent a hydrodynamic contribution that recovers "Einstein viscosity" at low z (dilute suspensions).Boyer et al. (2011) found out that all the data of n obtained for different suspensions collapse on a single n curve presented in Eq. ( 11) exhibiting that the friction coefficient n, i.e., the ratio eff x v , tends to a finite value 1 n at vanishing I v and increases with increasing I v .The existence of a finite n, a finite eff x v , at vanishing I v implies that the shear stress x of the dense suspensions has a "threshold" below which the shear stress cannot drag/haul the suspensions.Such a "yield condition" was also discovered in the case of dry dense granular flows and referred to as viscoplastic friction law (Jop et al. 2006).As for the volumefraction law which Boyer et al. (2011) proposed shown in Eqs. ( 9) and ( 12) above, the solids volume fraction z is a decreasing function of I v because when the shear strain rate c o increases (I v increases), sheared spherical solids will collide among each other to dilate (the porous volume expands).Consequently, z decreases.On the contrary, when c o decreases (I v decreases), initially "loose" solids will be compressed by the normal stress eff v to become more compact.In this case, z increases.As shown in Eq. ( 12), z increases to the maximum value m z at vanishing I v .By comparison, given a fixed eff v , the asymptotic behaviour of z close to o in the case of dense suspensions and is m \ z z c o for dry dense case.In our study here, we do not have laboratory data to describe the formulae of n and z for submarine landslides because the constitutive relations or material properties for submarine landslides are often inaccessible or imprecisely known.We assume that Eqs. ( 8) -( 12) obtained from the laboratory can be applied to describe submarine landslides.Accordingly, the shear stress x in Eq. ( 7) can be written as ( , ) ( ) ( ) 2.5 Therefore, Eq. ( 7) is re-written as where n(I v ) and I v are defined in Eqs. ( 15) and ( 16), respectively, and the shear strain rate c o is approximated by u/h.Note that in Eqs. ( 15) and ( 17), the friction coefficient n is time dependent because the excess pore pressure "p" in the definition of I v in Eq. ( 16) is increasing with time due to the basal frictional heating.This "increasing p" is discussed in section 3.2 below.Now, according to Eqs. ( 12) and ( 16), the rate of change of porosity in Eq. ( 6) can be expressed as where ( , ) cos I n p Substituting Eqs. ( 18) into (6), one can rewrite the governing equation for p as follows where 1 ( ) The first term on the RHS of Eq. ( 21) is governed by Eq. ( 17).
As shown in Eq. ( 17), the friction work xc o increases the temperature in the basal shear zone, which we call as "frictional heating".As a result, the volume of the pore fluid in the basal shear zone will expand in response to this heating, leading to an increase of excess pore pressure there.This "pore pressure increasing" phenomenon is called as "thermopressurization" (TP) below.Moreover, the thermopressurization coefficient Λ represents the amount of excess pore pressure generated per unit increase of temperature.The larger the Λ, the more excess pore pressure will be generated at a given rate of increase of temperature.The second term on the RHS of Eq. ( 21) is a "dilatancy" term.As one can see in Eq. ( 12), the experiment of Boyer et al. (2011) showed that when a dense suspension is sheared (c o increases and I v goes up), the sheared particles will collide among each other and therefore the porous volume will expand (the porosity increases and z decreases), which we call as "shear dilatancy" (SD).In this case, the water pressure within the expanding pore volume will decrease (Iverson et al. 2000).That is why a "negative sign" appears in front of the dilatancy coefficient }.
In addition to the TP and SD, it was discovered that a rapid sediment deposition can be a main cause of excess pore pressure.For example, we know that enormous numbers of very fine soil grains can be delivered by rivers or oceanic currents to rapidly deposit on the upper continental slope.Because such fine soil grains distinctively have high porous water contents (Hampton et al. 1996), if the deposition rate exceeds the expulsion rate of the porous water (i.e., the pore fluid does not escape fast enough during the rapid deposition), the accumulated weight of the deposited soil grains can consolidate/compress the porous water causing an increase of the pore water pressure to support a portion of the sediment overburden.Therefore, we assume that the total excess pore pressure (p) can be decomposed as where p s is the amount of excess pore pressure induced by sediment loading/deposition and pl is the rest part of p generated by the TP and SD (excluding sedimentation).However, the time and spatial scales of p s are, respectively, much longer and larger than those of pl .For example, p s often increases extremely slow in the time scales of months, years, or longer at regular sedimentation rates (normally several millimeters per year).Therefore, within tens of seconds, the rate of change of p s can be neglected (i.e., p s is viewed as a quasi-static term), , t p t 0 for time whithin tens of seconds Particularly, in this study, we only focus on the landslide initiation from t = 0 to tens of seconds, because (1) our objective is to identify whether a submarine landslide, once seismically triggered at t = 0, will continuously move downslope at t > 0 or slow down to a complete stop and because (2) "tens of seconds" is long enough to judge the result.Within this short duration, t = 0 -tens of seconds, the sedimentation-generated pore pressure, p s , can be viewed as a quasi-static term as described in Eq. ( 27).On the other hand, since the areas of sediment basins in the ocean are often very large (hundreds of km 2 ) and the sediment deposits in the basins are usually very deep (hundreds of meters), the spatial variations of p s in such sediment basins often have very long wave lengths.Therefore, within a basal shear zone the thickness of which is decimeters or smaller, the vertical gradient of p s can be neglected, Subsequently, by substituting Eqs. ( 26) -( 28) into ( 21), we obtain for 0 ≤ z ≤ h and 0 ≤ t ≤ tens of secs.The parameters Λ, }, A hy , and c o are defined in Eqs. ( 22) -( 25).Note that the superposition in Eq. ( 26) must be substituted into Eqs.( 14), ( 16), ( 17), and ( 19) to replace the total excess pore pressure "p".
In this article, we employ a finite-difference scheme (cf.Goren and Aharonov 2009) to approximate the secondderivative term on the RHS of Eq. ( 29) as The concept of Eq. ( 30) was also used in the study of fault slip (Segall and Rice 2006).In Eq. ( 30), l A t d h y d = is a diffusion length scale over which the excess pore pressure pl is relaxed, where t d is called "relaxation time".As suggested by Goren and Aharonov (2009) and a reference therein, the order of t d is about several seconds.The 'negative sign' in Eq. ( 30) represents that the excess pore pressure pl diffuses out of the shear zone.By substituting Eqs. ( 30) into ( 29) and setting z = 0, one can obtain the governing equation for ( , ) where ( , ) dT t dt 0 can be derived as follows by fixing z at z = 0 in Eq. ( 17) and using Eq. ( 26): In Eq. ( 32), we define the notation So, Eqs. ( 31) and ( 32) can be combined as The friction coefficient n in Eq. ( 34) will be extended to be a function of I v and temperature later, where I v , the viscous number, is re-defined as follows by setting z = 0 in Eq. ( 16) and using Eqs.( 26) and ( 33), ( , ) ( , ) cos c os The purpose of deriving Eq. ( 34) is to help us solve the landslide velocity u(t) in Eq. ( 1).Now, substituting Eqs. ( 26) and ( 33) into (1) to replace the total excess pore pressure p, one can re-write Eq. ( 1) as

Thermal Softening
Thermal softening behavior was discovered in laboratories (Hicher 1974;Despax 1976;Modaressi and Laloui 1997).For some kinds of clays, their friction coefficient decreases with the increase of temperature (see Fig. 1 in Vardoulakis 2002).Such property is called "thermal softening".Here, we assume the clay sediment in the shear zone has this "softening" property.By following the 'separationof-variables' expression in Eq. ( 24) of Vardoulakis (2002), Eq. ( 17) (Veveakis et al. 2007), and Eq. ( 20) (Veveakis et al. 2010), we assume that the friction coefficient n in our Eq.( 34) is of this form where Iv n is the same as Eq. ( 15) representing a "ratestrengthening" part and T n stands for the thermal softening part e ( ) The exponent M in Eq. ( 38) is a positive constant representing the decreasing rate of n per unit increase of temperature from a reference value T 0 .Here, we choose the value of M = 0.0093 (°C -1 ) suggested by Vardoulakis (2002) and Veveakis et al. (2007).As a matter of fact, how fast the friction coefficient decreases with temperature depends upon the clays (cf.Fig. 1 in Vardoulakis 2002).To summarize, in this study, we combine Eqs. ( 32), (34), and (36) as a system of ordinary differential equations (ODE) to study the occurrence of submarine landslides.The Gear-form fourthorder Runge-Kutta numerical scheme (Gear 1971) was used to solve this ODE system in a 'time-discrete step-by-step forward process'.By using this numerical method, at each time step t = t j , j = 1, 2, 3, …., the solutions of u(t j ), ( , ) p t 0 j l , and T(0, t j ) are obtained simultaneously (notice that t o = 0).

RESULTS AND DISCUSSIONS
As described in subsection 2.2, we assume that a large parallelepiped block of very fine soil grains overlies a gentle continental slope and suddenly at t = 0 the entire sediment deposit is seismically triggered to slide down the slope with an initial velocity.Our focus is to study whether this sediment block will keep moving over the gentle slope or slow down to a complete stop under the effects of the TPW, shear dilatancy, and thermal softening mechanisms.We study this problem by applying the fourth-order Runge-Kutta method to solve the system of Eqs. ( 32), (34), and (36).Before discussing our results, we introduce the parameter values used in our calculations.Unless otherwise specified, the parameter values defined in subsection 3.1 are fixed in this article.

Parameter Values
First, we introduce the range of p s (0, t).As mentioned in Eq. ( 26), p s denotes the pore pressure generated by sediment deposition, which can decrease the effective normal stress of the deposited sediments so the basal shear strength and the stability of the sediments could be reduced.In our study here, given a parallelepiped block of fine sediment grains that statically overlies a gentle slope as just described above, the larger the value of p s at the base (z = 0) of this sediment block, the smaller the basal shear strength of this sediment deposit.This large sediment block will be more unstable.That is to say, the entire sediment block is more unstable at smaller Qr [Notice that a larger p s (0, t) corresponds to a smaller Qr as defined in Eq. ( 33)].As one can see, the parameter Qr appears in the "effective normal stress" term on the RHS of Eqs. ( 32), (34), and (36).The lower Qr, the smaller the effective normal stress.The basal frictional resistance will decrease.Therefore, the parallelepiped sediment block can be more possibly triggered by smaller-magnitude earthquakes.Now, we discuss the ranges of p s (0, t) and Qr specified in this study.To date, it has been discovered that the magnitude of p s depends on sedimentation rates, soil permeability, and over-consolidated thickness (Gordon and Flemings 1998;Dugan andFlemings 2000, 2002;Sultan et al. 2004a;Flemings et al. 2008;Stigall and Dugan 2010).For example, as shown in the drilling measurements by Flemings et al. (2008) in the deepwater Gulf of Mexico, offshore Louisiana, the value of p s (called as "overpressure" in Flemings et al. 2008) increased with the vertically downward depth below the seafloor, e.g., increasing to 1 -2 MPa when the depth was increased to 200 -300 mbsf, and then the p s value changed insignificantly at depths > 300 mbsf.But the investigation of p s by Dugan and Flemings (2000) on the New Jersey continental slope showed a slightly different result.They found out that the value of p s , called as overpressure and denoted by p * in Dugan and Flemings (2000), kept an increasing tendency along the entire vertically downward depth; p s did not stop increasing even though the depth was increased downward to > 600 mbsf.As shown in Dugan and Flemings (2000, cf. Fig. 2 therein), the measured value of p s at depth 400 mbsf is about 2 -3 MPa caused by a sedimentation rate < 1 mm yr -1 .Accordingly, we believe that in a sedimentary basin with a much larger deposition rate (>> 1 mm yr -1 ) and a lower permeability (e.g., 10 -18 -10 -19 m 2 ), the magnitude of p s at a deeper depth below the seafloor (e.g., 400 mbsf) very possibly can exceed 3 MPa.Therefore, in our study here, given a large parallelepiped sediment deposit with a height D = 400 m ( t D gDcosi ≈ 3.997 MPa) overlying a slope of i = 1°, we assume that the sediment deposit was accumulated by a very rapid deposition (>> 1 mm yr -1 ) so that the magnitude of p s at the base (z = 0) of this deposit can be as large as 3.6 MPa, exceeding 90 % of the hydrostatic effective normal stress cos gD t i D . In our calculations below, we specify the range of p s (0, t) between zero and .c os gD 0 9 t i D , namely, Qr varying between one and 0.1.
In Eq. ( 1), we set f t = 1000 Kg m -3 , s t = 2700 Kg m -3 , L = 1000 m, W = 1000 m, D = 400 m, i = 1°, and C K = 0.002.The drag coefficient C D for the parallelepiped block is calculated by using the formula 1.95 -0.77D/W (De Blasio 2011) and the added mass coefficient add c calculated by (π/4 ~ 1.2) D/L (De Blasio 2011).Here, we take a lower estimation: add c ≈ 0.8 D/L.In Eqs. ( 22), ( 23), ( 25), and (34), we follow Table 1 (Goren and Aharonov 2009) to set f a = 10 -3 °C-1 , f b = 3 × 10 -9 Pa -1 , and f h = 10 -3 Pa s.Moreover, we assume that the permeability l in Eq. ( 25) ranges between 10 -17 to 10 -19 m 2 (cf.Ikari et al. 2009;Platt et al. 2014).The values of t in Eq. ( 36) and c t in Eq. ( 34) are determined in the following way.First, we assume a value of n (e.g., 0.416) and substitute it into the definitions of t and c t ( , c f = 4187 J Kg -1 °C-1 , c s = 1000 J Kg -1 °C-1 ) to get a "temporary value" of t and c t .Subsequently, we use the governing equations to calculate u(t) and ( , ) p t 0 l for t > 0 based on these temporary values of t and c t , and then substitute the solutions of u(t) and ( , ) for each time step t j , j = 1, 2, 3, ..., into Eq.( 35) to obtain the values of I v (t) at all these time steps t j , j = 1, 2, ..., etc. Next, we use Eq. ( 12) to calculate the variation of n(I v ) corresponding to this range of I v (t).If the value of n which we assume at the beginning (= 0.416) is within this range of n(I v ), we accept the temporary values of t and c t and the solutions of u(t) and ( , ) p t 0 l obtained based on these t and c t values.Otherwise, we choose a different value for n and re-do the same calculation in a try-and-error sense.
In this article, we choose two cases of h (the shear zone thickness), 0.001 and 0.1 m, to study the impacts of frictional heating on the initiation of submarine landslides.As shown in the studies of large terrestrial landslides, the value of h was assumed as a constant, e.g., 0.1 m (Goren and Aharonov 2009), 0.161 m (Veveakis et al. 2007), 0.02 m (Goren and Aharonov 2007), or some value between 0.0015 and 0.1 m (Vardoulakis 2002).Here, we compare those two cases of h (0.001 and 0.1 m) with the diffusion length scale, l d , to decide our simulation scheme.The order of l d can be evaluated by using the definition of l A t d h y d = shown near Eq. ( 30).We set t d to be of 2.35 and 20 s.The first one, 2.35 s, was suggested by Goren and Aharonov (2009) and a reference therein, whereas the second one, 20 s, selected in the following way.In this study, we calculated landslide evolutions only for time t = 0 -tens of seconds, which is long enough to judge whether a landslide will slow down to a complete stop or keep moving down a gentle slope.Therefore, we selected 20 s (one order larger than 2.35) as the second choice of t d .According to these two values of t d (2.35 and 20 s), we obtained that the order of l d is slightly less than 0.01 m.By comparing l d with those two h values (0.001 and 0.1 m), first, in subsection 3.2, where h is fixed at 0.001 m, in this case, because the value of l d much larger than h, implying that the pore pressure generated by the frictional heating diffuses out of the basal shear zone very much, the last term on the RHS of Eq. ( 34) is significant and should not neglected.On the contrary, in subsection 3.3, where h = 0.1 m, the value of l d much smaller than h, the last term on the RHS of Eq. ( 34) is insignificant and can be neglected.

Results for h = 0.001 m
Assume that a large parallelepiped block of fine sediment overlying a continental slope (i = 1°) was shaken by an earthquake to slide down the slope at an initial velocity, u(0) ≠ 0. The shear zone thickness is 0.001 m.For convenience, we use the abbreviation "IV" to stand for u(0).We set the initial conditions of ( , ) p t 0 l and T(0, t) to be zero and 12°C, respectively.For simplicity, the initial value of T(0, t) is chosen as the reference temperature T o in Eq. ( 38).The results for h = 0.001 m are exhibited in Figs. 1 -3.First, in Figs.1a and b, we show that given an IV, the sediment block with a larger (or smaller) Qr is more stable (or unstable).For instance, in Fig. 1a, given IV = 0.1 m s -1 , one can see that the sediment block keeps moving over the gentle slope for Qr = 0.1 (red line) but slows down to a complete stop for Qr = 0.5 (green) or larger.But if IV is increased to a larger value, e.g., 0.3 m s -1 as shown in Fig. 1b 1) and ( 36).As shown here, given an IV, the smaller the Qr (or given a Qr, the larger the IV), the more possible it is for landslides to continuously move down gentle slopes.
0.5, which was originally stable for IV = 0.1 m s -1 , becomes unstable and continuously slides down the gentle slope.Such a change of "slope stability" implies that the larger the earthquake, the more easily the block moves down the gentle slope (see Fig. 1c).Now, we select the case of IV = 0.3 m s -1 and Qr = 0.4 (the red line in Fig. 1b) to discuss the impacts of the TPW and thermal softening mechanisms on the initiation of submarine landslides.As shown in Fig. 1b, at initial, the slide velocity of the sediment block continuously decreases, but about t = 0.45 s, the block starts to accelerate moving faster and faster over the gentle slope.This is because combining the TPW and thermal softening mechanisms can quickly decrease the basal frictional resistance of the block to a degree less than the "slope-parallel" component of the effective weight of the block, so that this sediment block, before it stops, can change to acceleration continuously moving down the gentle slope.As this block slides downslope, the basal frictional stress does a mechanical work to increase the internal thermal energy of the water-sediment mixture in the basal shear zone.Consequently, as we present in Fig. 2a, the basal temperature, T(0, t), increases with time.This is the frictional heating mentioned above.Furthermore, as shown in Fig. 2b, the friction coefficient n decreases with time due to the thermal softening effect.Next, we exhibit the impacts of thermopressurization, shear dilatancy, and diffusion on the generation of the excess pore pressure ( , ) p t 0 l . First of all, we introduce three abbreviations.The first one, TP, represents the first term on the RHS of Eq. ( 34 The other two, SD and DIF, stand for the second and third terms on the RHS of Eq. ( 34), respectively.In other words, TP, SD, and DIF represent the thermopressurization, shear dilatancy, and diffusion contributions to ( , ) p t 0 l , respectively.In Fig. 3a and the inset, we use red, green, and blue lines to represent the temporal variations of TP, SD, and DIF, respectively.As shown, the value of TP decreases rapidly from its maximum (= 14.1 MPa s -1 ) occurring at t = 0 down to an equilibrium value (= 0.59 MPa s -1 ) within about one second.Such rapid decrease of TP shows that the excess pore pressure ( , ) p t 0 l increases very fast.Moreover, the temporal evolutions of SD (green) and DIF (blue) are similar to the trend of TP (red), the details exhibited in Fig. 3b.As shown there, the SD (or DIF) line decreases rapidly from the maximum value to the equilibrium one, about -0.01 MPa s -1 (or -0.58 MPa s -1 ), within only one second.Obviously, during the period when ( , ) p t 0 l increases very fast, about t < 0.5 s, the contributions of SD and DIF to the generation of ( , ) p t 0 l are negligible compared with that of TP; the contribution of TP dominates the generation of ( , ) p t 0 l during this period.Moreover, about t = 0.45 s, the sediment block changes to acceleration and then keeps moving down the gentle slope as we explain in Fig. 1b.This is the evidence showing that the TP term is a main cause of the decreasing of the basal frictional resistance.Because of the TP contribution, the sediment block can overcome the basal frictional resistance to continuously slide down the gentle slope.During t > 1 s, as shown in Fig. 3b, the DIF term nearly balances the TP one.Note that the value of SD changes from positive to negative at t = 0.45 s because at this moment the downslope velocity of the block changes from deceleration to acceleration.
Lastly, there are two minor things that need explanations.First one, the magnitude of the water resistance, the last term on the RHS of Eqs. ( 1) and ( 36), is about 1/100 -1/1000 of the magnitude of the first term during t ≤ 25 s, so this "water resistance" can be neglected during t ≤ tens of seconds.Secondly, the calculated values of the porosity (n) vary between 0.4155 and 0.4165 for t ≤ 25 s, this n(I v ) range covering the value of n which we assume at the beginning (= 0.416).

Results for h = 0.1 m
The results for h = 0.1 m are exhibited in Fig. 4. Here, we let all the parameter values, except h and IV, be the same as those specified in Fig. 3.In Fig. 4, the shear zone thickness h is assumed to be 0.1 m, which is 100 times as large as the one (= 0.001 m) defined in Fig. 3. Our calculated results in Fig. 4 exhibit that the landslide continuously slides down the gentle slope at IV ≥ 1.52 m s -1 .But if IV < 1.52 m s -1 , e.g., 0.3 m s -1 , the landslide will decelerate to stop.On the contrary, the landslide discussed in Fig. 3, with the same IV (= 0.3 m s -1 ) but smaller h (= 0.001 m), continuously slides down the same gentle slope.In theory, the larger the thickness "h", the smaller the TP as shown in Eq. ( 39).Meanwhile, the smaller TP, the lower the mobility of landslides as explained in subsection 3.2.Therefore, the case of landslide discussed in Fig. 4 with h = 0.1 m has lower mobility than the one in Fig. 3 with h = 0.001 m and thus more easily decelerates to stop.

CONCLUSIONS
Up to date, it is still an open question why submarine landslides can travel far on gentle slopes (Talling et al. 2007).In our here, we numerically show that the frictional heating in the basal shear zone can "lubricate" the movement of submarine landslides.The larger the heating, the higher the slide mobility of submarine landslides.With the help of this "frictional heating lubrication", submarine landslides could slide farther even on very gentle slopes.

Fig. 3 .
Fig. 3. Temporal evolutions of TP (red), SD (green), and DIF (blue).As shown in part (a) and (b), at t < 0.5 s, the effects of SD and DIF are negligible; TP dominates the generation of ( , ) p t 0 l during this period.About t > 1 s, DIF nearly balances TP; their magnitudes are about the same ≈ 0.5 Mpa s -1 .Note that all the parameter values and initial conditions specified in Figs. 2 and 3 are the same.

Fig. 4 .
Fig.4.Landslide velocity vs. time for h = 0.1 m, IV = 2.00 (red), 1.52 (green), and 1.00 (blue) m s -1 .All other parameter values are the same as those specified in Figs.2 and 3. Our calculated results illustrate that IV = 1.52 m s -1 is a critical value, above (or below) which the sediment block continuously moves downslope (or decelerates to a complete stop).

Table 1 .
Variables, parameters, and abbreviations used in section 2.