A general analytical solution for the evolution of cliffs accounting for strength degradation, seismic action, formation of tension cracks and seepage

Article history: Received 4 June 2016 Received in revised form 16 September 2016 Accepted 18 September 2016 Available online 20 September 2016 The evolution of natural slopes over time is ruled by several concurrent physical phenomena, namely the strength of its component geomaterials and their weakening over time due to weathering processes, the occurrence of seismic events, seepage and the formation of tension cracks. The paper presents analytical solutions obtained considering a succession of discrete failure events (landslides) progressively altering the slope morphology over time. The model, derived in the framework of limit analysis assuming plane strain conditions, provides a tool for the assessment of whether manufacts and/or infrastructures located on a slope subject to various natural degradation phenomena will be affected by the occurrence of failures. Unlike current empirical and semi-empirical models of slope evolution, the analytical solution that is here presented is derived by applying principles of soil and rock mechanics, therefore it is of general validity, so that no ad–hoc calibration against past observations of the evolving slope is needed. This analytical technique only requires knowledge of the (geotechnical) parameters characterising the geomaterials comprising the slope of interest, namely angle of shearing resistance, φ, cohesion, c, tensile strength, unit weight, together with knowledge of the relevant seepage scenarios, strength degradation processes, and seismic events likely to occur. Results show that that earthquake loading and seepage can substantially decrease slope stability, increase the volume of material sliding away during each landslide event and alter the evolution of the slope over time, while tensile strength exhibits a less strong influence especially as strength degradation progresses. Dimensionless ready-to-use charts are provided for the benefit of practitioners. © 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Morphological evolution of cliffs is a traditional subject in engineering geology and geomorphology (Carson and Kirkby, 1972;Hutchinson, 1969Hutchinson, , 1970Hutchinson, , 2001Selby, 1982;Parsons, 2002). Modelling the progressive retreat of cliffs has recently received considerable attention by the engineering community due to increasing coastal erosive processes caused by climate change and increased environmental awareness at national and European level (Bray and Hooke, 1997). Also the insurance industry needs reliable models for the predictions of the amount of cliff retreat over time for residential buildings located in exposed areas whereas local authorities and decision makers need to know the level of risk faced by residential buildings and public infrastructure (e.g. coastal roads, pedestrian footpaths, car parks, etc.).
Cliff evolution is caused and/or accelerated by several physical agents (Arkin and Michaeli, 1985;de Lange and Moon, 2005;Briaud, 2008;Collins and Sitar, 2011). Key drivers of slope instability are seismic action (Chen and Liu, 1990;Ling and Leshchinsky, 1995;Loukidis et al., 2003;Wasowski et al., 2011;Rathje and Antonakos, 2011;Yang and Chi, 2014;Tsai and Chien, 2016), rainfall and climatic variations (Leroueil, 2001;Frayssines and Hantz, 2006;Take and Bolton, 2011;Conte and Troncone, 2012;Springman et al., 2013), weathering (Yokota and Iwamatsu, 2000;Hachinohe et al., 2000), crack formation (Baker, 1981;Hales and Roering, 2007) and wave action for seacliffs (Benumof et al., 2000;de Lange and Moon, 2005). In this paper we shall consider all these actions but the last one. The first stability analyses of slopes were based on limit equilibrium methods (Espinoza et al., 1994). However, they present several shortcomings (Duncan, 1996) one of which is the fact that the solution is neither a lower nor an upper bound of the true collapse load. The use of numerical methods such as finite element method (Zheng et al., 2005;Huang and Jia, 2009;Potts and Zdravković, 2001) and material point method (Yerro et al., 2015) to provide approximate solutions to the slope stability problem is also increased in the latest decades. The Discrete Element Method has also been employed (Utili Engineering Geology 219 (2017)

Contents lists available at ScienceDirect
Engineering Geology j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / e n g g e o Notation A area of the sliding wedge c cohesion d 1 , d 2 lengths as defined in Fig. 3 f 1 ,f 2 , etc. mathematical function for the external work rate done by the gravity f d function for the dissipated energy f 1v , f 2v , f 3v , etc. mathematical function for the external work rate done by the gravity, f 1h , f 2h , f 3h , etc. mathematical function for the external work rate done by the horizontal seismic inertia F PSh horizontal pseudo-static force F PSv vertical pseudo-static force f t unconfined tensile strength measured in laboratory experiments f t M−C unconfined tensile strength predicted by the Mohr-Coulomb criterion G weight of the failing wedge g gravitational acceleration H slope height K h horizontal seismic acceleration coefficient K v vertical seismic acceleration coefficient L, L 1 , L 2 , l 1 , l 2 lengths as defined in Fig. 1 m 1 , m 2 masses of the sliding wedges N s stability factor r generic radius of curvature of the logarithmic spiral r x , r y , r z radii of the spiral at the angles x, y and z r u pore pressure coefficient _ u displacement rate (vector) Ẇ ext rate of external work Ẇ d rate of internally dissipated energy Ẇ cf rate of energy dissipated by crack formation Ẇ 1, Ẇ 2, etc. rates external work x 1 , y 1 , z 1 , x 2 , y 2 , z 2 , μ 1 , μ 2 angles as defined in Fig. 1 β inclination of the slope front γ soil unit weight δ tension crack depth θ generic angle of the logarithmic spiral θ 1 , θr u 1 , θ 2 , θr u 2 angles as defined in Fig. 3 ω angular velocity ϕ angle of shearing resistance and Crosta, 2011b; Camones et al., 2013;Utili et al., 2015) and also for 3D analyses of the stability of rock slopes (Boon et al., 2014) which has been made possible by computational advances in the DEM contact detection algorithms to deal with polyhedral blocks (Boon et al., 2012(Boon et al., , 2013 and new algorithms for rock slicing (Boon et al., 2015). However, in cases of rather uniform slopes subject to a few cracks, such a numerical approach is not justified. This paper aims to develop an analytical solution describing the morphological evolution of natural cliffs with limited tensile strength subject to progressive retreat induced by ground strength degradation, for instance weathering induced, for static and seismic scenarios and various hydrological conditions. The limit analysis upper bound method (Chen, 1975) and the pseudo-static approach (Terzaghi, 1943) are used to derive the analytical law describing the evolution of homogeneous c, ϕ slopes subject to strength degradation, seismic action, formation of tension cracks and various seepage conditions. Perhaps the greatest limitation of the analytical solution presented resides in assuming the slopes homogeneous. In fact natural heterogeneity, layering, different material properties of the layers etc. tend to produce geomorphic features unique to each particular slope that are not captured by the solution here presented. However, the other side of the coin is the fact that such a strong assumption allows for the derivation of a comprehensive analytical solution that can be used to achieve a first rough estimation of the past or future evolution of a slope knowing a limited amount of information and to explore the relative influence on slope evolution of the various physical phenomenon considered in general terms.
The ground strength is here characterised by the Mohr-Coulomb failure criterion therefore only three parameters are needed to describe its properties (unit weight, angle of shearing resistance and cohesion). The limit analysis upper bound method was applied to determine each discrete landslide event occurring over time for successive destabilization and complete removal of the failed mass after each event (Utili and Crosta, 2011a). The material accumulated at the slope toe cannot be taken into account in our model since limit analysis is not able to give any information about the final geometry of the debris accumulated after each landslide. Therefore, it is assumed that the debris accumulating at the slope toe is removed by atmospheric agents and erosion (e.g. fluvial or marine), before a new landslide develops. This condition is known in the literature as a strong erosion condition and is typical of weathering-limited processes (Hutchinson, 1973).

Limit analysis model
The failure mechanisms assumed in our analysis are 2D single wedge rigid rotational mechanisms (see Fig. 1). The failing wedge EBCD rotates rigidly around point P 1 with the ground lying on the right of the log-spiral CD and the vertical crack BC remaining at rest. The equation of the logspiral CD is: with r being the distance of a generic point of the spiral to its centre (see Fig. 1).
For sake of simplicity, the calculations presented below refer to the case of a horizontal slope crest. Since the limit analysis method does not provide information about the final geometry of the debris accumulated after each landslide, it is therefore assumed that the debris is entirely removed by atmospheric agents and/or erosion before the occurrence of the next slope failure.
The following geometrical relationships will be employed in the derivation of the semi-analytical solution: L 1 ¼ r x1 cos x 1 − cos y 1 exp tanϕ y 1 −x 1 ð Þ ½ ð7Þ L 2 ¼ r x 2 cos x 2 − cos y 2 exp tanϕ y 2 −x 2 ð Þ ½ ð9Þ where r x, r y and r z are the radii of the spiral at the angles x, y and z respectively, L, L 1 , l 1 L 2 , l 2 are the horizontal lengths (see Fig. 1) and H is the height of the slope.

First failure (landslide)
The calculations for the first failure mechanism are here illustrated. The upper bound is derived by imposing energy balance for the failing wedge EBCD: where Ẇ d is the rate of energy dissipated along the log-spiral surface (shear failure), Ẇ cf is the rate of energy dissipated by crack formation and Ẇ ext the rate of external work done by the failing wedge. The calculation of Ẇ d accounting for the energy dissipated along the log-spiral segment CD is reported in (Utili, 2013). Cracks may develop from the slope face (DE) and/or from the upper part of the slope (EF) (see Fig. 1). Cracks or fissures can be the result of a variety of phenomena, for instance exceedance of the ground tensile strength (Baker, 1981), the occurrence of differential settlements (Vanicek and Vanicek, 2008), desiccation (Konrad and Ayad, 1997;Dyer et al., 2009;Péron et al., 2009) and freezing (Hales and Roering, 2007). Cracks may cause a significant decrease in the stability of a slope both in static (Michalowski, 2013;Utili, 2013;Gao et al., 2015;Utili, 2015) and seismic conditions (Utili and Abd, 2016).
Here cracks are treated as no-tension non-cohesive perfectly smooth (no friction) interfaces, therefore the angle η between the velocity vector and the crack surface (see Fig. 1) is 0°b η b 180°. Following (Michalowski, 2013), two types of cracks can be considered: cracks existing in the slope before the formation of any failure mechanism, called pre-existing cracks by Michalowski, and cracks forming as part of the failure process due to the exceedance of the ground tensile strength, here called tension cracks, that take place contemporaneously to the formation of localized deformations leading to the failure of the slope. A pre-existing crack may have been formed in the past by tensile stresses that are no longer acting on the slope, for instance due to tectonic movements. Tension cracks instead indicate cracks that are generated by tensile stresses exceeding the current ground tensile strength leading to the formation of a failure mechanism (Terzaghi, 1943;Baker, 1981). Although the presence of pre-existing cracks can be easily accounted for in limit analysis (Michalowski, 2013;Utili, 2013), here only tension cracks are considered. This is because the evolution of a slope subject to deterioration/weathering is ruled by the deterioration of the ground tensile strength, leading to the onset of tension cracks whereas pre-existing cracks generated in the course of the geological history of the slope formation are likely to affect the first failure mechanism only. To calculate the energy dissipated by the formation of a crack with limit analysis, Michalowski (2013) has considered limiting the Mohr Coulomb linear envelope by the stress circle of an unconfined uniaxial tensile strength test with the circle being tangent to the linear envelope (see Fig. 2). This failure criterion is realistically non-linear in the tension zone and on the other hand lends itself to limit analysis calculations. The energy expended for the formation of a tension crack, Ẇ cf turns out to be (Michalowski, 2013): with μ 1 being the angle made by the segment P 1 B with the horizontal (see Fig. 1), f c M−C being the Mohr-Coulomb unconfined compressive strength of the ground and f t the unconfined tensile strength as measured from laboratory experiments (see Fig. 2). It is convenient to introduce a dimensionless coefficient, t, defined as the ratio of the unconfined tensile strength measured in laboratory experiments, f t , over the full unconfined tensile strength predicted by the Mohr-Coulomb criterion, f t M−C (see Fig. 2a): It is straightforward to observe that 0 b t b 1. Also both f c M−C and f t M−C are uniquely related to c and ϕ: Three different conditions controlling crack formation are tackled here: i. tension cut-off (t = 0); ii. limited tensile strength, e.g. t = 0.2 and t = 0.5; Fig. 4. Potential failure mechanisms for the second failure, relative to different mechanisms considered (for different critical heights h i ). iii. full tensile strength, i.e. the possibility of tensile failure of the ground is neglected The rate of external work for the sliding wedge EBCD, Ẇ ext , is made by two contributions: with Ẇ γ representing the external work done by the weight of the wedge and Ẇ w representing the work done by the water pressure. Ẇ γ is calculated as the work of block EFD minus the work of block BFC. In turn the work of block EFD is calculated by algebraic summation of the work of blocks P 1 FD, P 1 FE and P 1 ED (Chen et al., 1969) and the work of block BFC is calculated by summation of the work of blocks P 1 FC, P 1 FB and P 1 BC (Utili and Nova, 2007;Utili and Crosta, 2011a;Utili, 2013). Note that here, in addition to the weight force, a horizontal pseudo-static force, F PSh = mK h g = γK h A, with g being the gravitational acceleration and m the mass of the wedge, and a vertical one, F PSv = mK v g = γK v A, are added to account for seismic action (Chen and Liu, 1990;Utili and Abd, 2016).   To calculate Ẇ w , a uniform value of the dimensionless coefficient r u is assumed over the whole slope with r u defined as (Bishop and Morgenstern, 1960): with u being the total pore water pressure in the considered point of the failure line, γ sat the ground bulk unit weight and h the depth of the point considered from the ground surface. The assumption of uniform r u is a strong one (Barnes, 2010), however this assumption is still commonly made in slope stability analysis (e.g. Michalowski, 2013). The work of water along the cracked part BC and the log-spiral part CD (Fig. 3) was calculated as a surface integral over the whole surface (Michalowski, 1995(Michalowski, , 2013: where [v] n is the normal component of the boundary velocity (Fig. 3) and tan θr Substituting Eqs. (12), (16) and (18) into Eq. (11), the final equation to calculate the stability factor, N = γΗ/c, is obtained: with λ = K v /K h (consistently with Fig. 1, the + sign indicates vertical downward acceleration, whereas thesign indicates vertical upward acceleration). The global minimum of g (x 1 , y 1 , z 1 , ϕ, β, K h , λ) over the three geometrical variables x 1 , y 1 , z 1 provides the least (best) upper bound on the stability factor for the case that has been considered. The static case is a particular case obtained setting K h = K v = 0. The calculations of the expressions for Ẇ ext for each block are provided in Appendix A.
(a) r u =0.25 (b) r u =0. 5  Unlike the case of intact slopes, failure mechanisms may in principle daylight on the slope face above the slope toe; however in (Utili and Abd, 2016) no potential mechanism passing above the slope toe turned out to be the most critical, therefore no potential failure mechanisms passing above the toe were considered in this analysis.

Second and successive failures (landslides)
The analytical expressions for the second failure also apply to every successive failure. After the region EBCD (Fig. 1) has slipped away and due to further ground degradation, at some point a second landslide will occur. The double logarithm spiral shaped area GIDCB will rigidly rotate around the center of rotation P 2 , as yet undefined, with the material below the logarithmic spiral ID and right of the vertical crack GI remaining at rest. The mechanism is now defined by six variables x 1 , z 1 , y 1 , x 2 , z 2 , y 2 , where x 1 , y 1 , z 1 are the angles defining the first log spiral failure line (i.e., the current slope  profile produced by the previous failure), that will be called 'old' landslide and will be denoted by the superscript o and x 2 , y 2, z 2 are the angles defining the second log spiral failure line (see Fig. 1), that will be called 'new' landslide and will be denoted by the superscript n.
Two contributions constitute the rate of the external work for the sliding wedge GIDCB; the work done by the weight of the wedge that slides away Ẇ γ and the work done by the water pressure Ẇ w . The weight of the wedge that slides away Ẇ γ is calculated this time as the work by block KJD minus the work of blocks KBCD and JIG. The work of block KJD is calculated as the summation of three contributions:  Values for associated crest retreat and sliding area for slopes with initial inclination β and angle of shearing resistance ϕ, for slopes with t = 1.0, t = 0.5, t = 0.2 and t = 0.  Table 2 Values for associated crest retreat and sliding area with initial inclination β and angle of shearing resistance ϕ, for intact slopes under seismic action or with the existence of water (t = 1). Ẇ 2 n and Ẇ 3 n being work rates done by blocks P 2 JD, P 2 JK and P 2 KD respectively. The rate of work of block KBCD is calculated as the summation of six o being the work rates of regions P 1 FD, P 1 FK, P 1 KD, P 1 FC, P 1 FB and P 1 BC respectively. Τhe work rate of JIG is calculated as Ẇ 4 n -Ẇ 5 n -Ẇ 6 n , with Ẇ 4 n , Ẇ 5 n and Ẇ 6 n being the work rates done by P 2 JI, P 2 JG and P 2 GI respectively. Therefore, Ẇ γ is calculated as: with n and o referring to the new and the old landslide respectively. The calculations of the expressions for Ẇ 1 n , Ẇ 2 n , Ẇ 3 n etc. and for Ẇ w are provided in Appendix A. Substituting Eqs. (12), (18) and (20) into Eq. (11) and dividing all terms by _ θ and r x 2 , and rearranging, the stability factor, N s = γH/c, is obtained: The global minimum of g (x 1 , y 1 , z 1 , x 2 , y 2 , z 2 , ϕ, β, K h , λ) over the three geometrical variables x 2 , y 2 , z 2 , provides the least (best) upper bound on the stability factor for the second and any successive landslide. Note that the second (and every successive) mechanism could pass through any point since the current slope profile is no longer straight, as presented in Fig. 4. Therefore failure mechanisms daylighting at any point of the slope profile left after the first failure has occurred have to be considered. To this end, the slope profile was divided into a discrete number of points (n) and each point has been assumed as the toe of a sub-slope whose height, h i , is smaller than the overall height H (see Fig. 4). The most critical mechanism among all the possible mechanisms has to be found. The critical cohesion value, c i , and angles, x i , y i and z i , associated with the critical failure mechanism, were determined for all n sub-slopes of different height, h i , with the parameter y i assuming a different value for each sub-slope analysed. The most critical failure mechanism among the n potential mechanisms is the one associated to the highest cohesion value. As long as a sufficiently large value of n is considered, n does not affect the obtained results.

Weathering induced cohesion decrease
By applying the procedure described, it is possible to determine as many failure mechanisms as needed to follow the slope evolution until full degradation of the ground strength has taken place. Here, the friction angle is assumed to remain constant in time with only cohesion decreasing over time, since experimental evidence from the weathering of rocks and cemented soils show that weathering causes mainly a decrease of cohesion and to a much lesser extent to the angle of shearing resistance (see Fig. 5). However, extension of the model to account for a decrease of both c and ϕ is straightforward. The case of both c and ϕ decreasing is reported in (Utili and Crosta, 2011a) for geomaterials with full tensile strength.
Computations for the evolution of slopes have been carried out using Matlab for a wide range of parameters (ϕ and initial slope inclination β), for intact slopes (t = 1), slopes with soil tensile strength limited to t = 0.2 and t = 0.5 and slopes with no tensile strength, t = 0 (tension cut-off) under different seismic and seepage scenarios.

Morphological evolution of slopes
In Fig. 6 the typical evolution undergone by a slope subject to weathering induced strength degradation is shown. The features of the morphological evolution exhibited by the slope are similar to what shown in (Utili and Crosta, 2011a) apart from the upper part of the failure mechanisms which is always vertical due to the presence of cracks. A number of successive failures of decreasing area and depth of mechanism occur until a deeper mechanism of much larger area takes place (see Fig. 6). To calculate the deeper mechanism, the slope profile, which is composed of several log-spiral pieces, was approximated by a straight line obtained as the linear envelope of the piecewise log-spiral profile (see Fig. 6). This approximation is acceptable since the region delimited by the log spiral and the vertical line is small in comparison with the size of the failure mechanisms.
In Fig. 7, the evolution of initially steep slopes, with very limited or completely absent tensile strength, is shown. It emerges that this type of slope exhibits a peculiar type of evolution: several successive small failure mechanisms made of thin slices occur (from vertical line FC to vertical line GI) until a much larger failure mechanism takes place (mechanism JMK). The occurrence of thin slice mechanisms, which in the limit case of an initially vertical slope (β = 90°) become infinitesimally thin slices, is described and physically explained in (Utili, 2013).

Parametric analysis of slope evolution for various hydrogeological scenarios
In Fig. 8 the effect of the seismic action is analysed for various values of tensile strength. For each failure, numbered in chronological order, the area of the sliding mass is plotted. It emerges that the influence of the soil tensile strength on the sliding area is small, but seismic action has the effect of amplifying these small differences. To investigate the influence of seepage instead, the critical height for a slope with an initial inclination of β = 60°and ϕ = 20°is plotted in Fig. 9 for 6 successive failures.
In Fig. 10, the step-like relationship between dimensionless normalised cohesion and crest retreat is plotted for 8 successive failures for various seismic scenarios (K h = 0, K h = 0.1 and K h = 0.2). It can be observed that if the first two failures are excluded, the values of critical cohesion and crest retreat lie on straight lines in agreement to what observed in (Utili and Crosta, 2011a) for static dry slopes not subject to crack formation. From the figure it also emerges that the tensile strength does not affect significantly either the cohesion or crest retreat. However this is less the case for slopes subject to substantial seismic action.
The evolution of the stability number (γH/c) and the dimensionless cliff retreat (L/H) as a function of the friction angle represented in the same graph as curves from the analysis of six successive failures, for initial slope inclinations of β = 60°, β = 70°and β = 80°are plotted in Fig. 11. Solid black lines are used for the stability number and dashed lines for the dimensionless crest retreat, while the corresponding number of failure is reported on each curve. It is evident that the influence of tensile strength and of the presence of cracks on the stability of the slopes is more significant in steep slopes. When seismic acceleration is taken into consideration the influence on both stability number and cliff retreat is significant (Fig. 12). It can be seen in Fig. 13 and as explained by (Michalowski, 2013), that the adverse influence of pore-water pressure on the stability number increases for larger values of ϕ.
Note that for some cases of slopes with low ϕ, after a number of failures, the model predicts a much larger landslide in terms of area (volume), but not in terms of crest retreat (Fig. 6) so this is more of interest for anthropic constructions at the toe of the slope. If the slope is a cliff by the oceanan environment where continuous wash away actions take placeit is reasonable to expect the erosion of the failed material before the next mechanism occurs. Instead, in case of less exposed cliffs (e.g. inland) we can expect that this deep mechanism will not take place due to the stabilising action of the weight of the debris and therefore it is reasonable to disregard it.
In Table 1 values of crest retreat and area of the failure mechanism are reported for various initial slope inclinations (β = 60°, 70°and 80°), angle of shearing resistance (ϕ = 20°, 30°and 40°) for the three scenarios of crack formation considered in the paper (tension cut-off, limited tensile strength and full tensile strength). The results of the parametric investigation for initially intact slopes subject to seismic action and various seepage scenarios are reported in Table 2.

Conclusions
Analytical solutions have been presented to investigate the effect of weathering induced strength degradation, seismic action, crack formation and various seepage scenarios on the geomorphological evolution of homogeneous slopes, employing the upper bound limit analysis method. Slopes are characterised by uniform cohesion, angle of shearing resistance, and finite tensile strength. The solutions were obtained considering a succession of discrete failure events (landslides) progressively altering the slope morphology over time and can be used to achieve a first rough estimation of the past or future evolution of a slope knowing a limited amount of information.
A comprehensive parametric analysis has been carried out to investigate the relative influence on slope evolution of the physical phenomena considered, i.e. weathering induced strength degradation, seismic action, crack formation and seepage. It emerges that strength degradation, seismic action and seepage exhibit a stronger influence on the morphologic evolution of slopes than the formation of cracks. High seismic acceleration and/or seepage forces cause larger slope failure mechanisms and therefore larger slope crest retreats. Tension cracks forming due to limited tensile strength influence the first slope failures more than the subsequent ones.