Investigation of the role of crown crack in cohesive soil slope and its effect on slope stability based on the extended finite element method

Tensile cracks in soil slopes, especially developing at the crown, have been increasingly recognized as the signal of slope metastability. In this paper, the role of crown cracks in natural soil slopes was investigated and its effect on stability was studied. A numerical modeling of slope and simulation of tensile behavior of soil, based on the Extended Finite Element Method (XFEM), was used. A numerical soil tensile test was applied to validate the use of XFEM on tensile behavior of soil before the simulation. Slope failure was simulated by using the Strength Reduction Method, which determines the potential slip surface of slope. The simulation results indicate forming of crown crack in natural soil slopes when the plastic zone starts penetrating. Therefore, it is reasonable to consider the crown crack as the signal of slope metastability. A sensitivity analysis shows the effect of cracking on slope stability if cracks are at the position of the tension zone. The stress variation analysis, from the surface slip deformation, reveals that the slope is at a state of compressive stress. When plastic zone starts penetrating, the upper part of slope generates tension zone, but the extent of tension zone is restricted by the slope failure. This suggests why tensile cracks are difficult to form and stretch in the deep part of the slope. The implementation of XFEM on slope stability analysis can be used for assessing the tensile strength of soil and forecasting the time of slope failure related disaster.


Introduction
Landslides are one of the most common geological natural hazards. Every year, many disastrous landslides occur all over the world, causing loss of property and life (Tang et al. 2009;Yin et al. 2016;Gianvito et al. 2018;Chen et al. 2018b;Zhang et al. 2019;). Landslide is a type of mass movement, which includes diverse ground movements such as fall, topple, slide, spread and flow (Hungr et al. 2014). Most of these ground movements indicate mass separation behavior, generating cracks prior to mass separation. Some landslides reveal the whole process of failure Ouyang et al. 2019;Chen et al. 2019), and a typical example would be the Baige landslide on the Jinsha river in China. Researchers applied remote sensing images and Interferometric Synthetic Aperture Radar (InSAR) to the Baige landslide, to analyze its historical deformation ; (Ouyang et al. 2019). They found that the swelling firstly occurred at the toe of the landside and, as time goes on, large-scale tensile cracks and drop heads formed at the crest of the slope. Finally, with continuous deformation, a penetrating slip surface generated, and a slide occurred. This case shows a common landslide failure process. Although, many events of landslides (Steiakakis et al. 2009;Chen et al. 2018a;Tang et al. 2019;Fan et al., 2019;Zhang et al. 2021a;Zhang et al. 2021b) and field investigations, carried out by the authors (Fig. 1) have proven that crown cracks are the signal of slope metastability during the long-term slow deformation process, there are still some issues to be investigated, such as: at what stage of slope deformation will the crown cracks form? How do these cracks form? What's the maximum length that cracks can be? Do these cracks affect the slope stability? Although laboratory tests have been carried out to answer some of these landslide issues (Tang et  , some data cannot be easily acquired in this way, which restricts further research on the matter. Why the depth of a crack cannot develop downward without limitation and other questions remain unanswered even though knowing the answers is important for comprehending the landslide behavior and prevention. Numerical simulation can help researchers obtain extra data, which are not available in laboratory tests. Currently, the Limit Equilibrium Method (LEM) and the Finite Element Method (FEM) are the most common methods to evaluate slope stability with cracks (Bishop 1973;Seed et al. 1990;Griffiths and Lane 1999;Qu et al. 2009;Bao et al. 2019;Lei et al. 2021). Additionally, Discrete Element Method (DEM) (Zhou et al. 2009a, b;Bao et al. 2020), Smoothed Particle Hydrodynamics (SPH) Ray et al., 2019) and Material Point Method (MPM) Conte et al. 2019) have been developed for analyzing slope stability as well as the post-failure movement. Using these new methods for slope stability analysis depends on the Strength Reduction Method (SRM) and they have advantages for the analysis of large deformation issues. However, DEM, SPH and MPM still have some limitations in stability analysis. In DEM, the material behavior depends on interactions between particles. It is challenging to microscopically measure some variables between particles and there is no well-substantiated theory to illustrate the relation between the magnitude of macroscopic and microscopic parameters. Thus, it is time-consuming to determine proper variables, especially in the SRP. The SPH and MPM are also time-consuming compared to LEM and FEM. The computation accuracy of SPH and MPM is lower than FEM within small deformation process prior to slope failure. Even though in theory these new methods can be used for slope stability analysis, they lack significant engineering validation compared to LEM and FEM, especially when involving cracks. As for the LEM, it is a classical method which has been used for the slope stability analysis involving cracks (Tang et al. 2019). However, there are many restrictions within LEM to analyze a problem with cracks as LEM is not rigorous and is limited because location and geometry of cracks need to be assumed in the model, and cannot be updated with calculations. Some researchers (Michalowski 2013;Utili 2013;Gao et al. 2015) develop Limit Analysis (LA) on stability of slopes with cracks, based on rigid mechanism and can overcome deficiency that the crack needs to be initialized in the LEM. LA can be used to explore cracks of unknown depth and unspecified location. Nonetheless, LA is mainly used to study slopes with simple geometry . FEM can subdue above limitations, which makes it a potential choice for slope stability analysis involving cracks.
The conventional FEM (CFEM) is difficult to use for implementing a highly nonlinear problem (Geng et al. 2018(Geng et al. , 2021. Simulating discontinuous elements, such as cracks, is also challenging due to meshing limitation. To overcome the restrictions of CFEM in the discontinuous analysis, some theories are proposed, including the advanced remeshing techniques (Areias et al. 2013(Areias et al. , 2015, the Numerical Manifold Method (NMM) (Shi 1991;Ma et al. 2009) and the extended FEM (XFEM) (Moës et al. 1999). Most of the methods modeling crack propagation by FEM heavily depend on the mesh alignment (Rabczuk and Ren 2017), while the XFEM can avoid the problem. In XFEM, special functions and element segmentation methods are used to fuse the solution of finite element approximation. The Level set Method (LSM) is used to show the geometry and extension process of the discontinuous interface. Differing the CFEM, XFEM no longer has rigorous requirements on the accuracy and repetition of the network, no specific restrictions on the crack front and growth path and has a high computational efficiency. Therefore, this method has been widely applied in fracture mechanics and engineering (Sanborn and Jean 2011;Wang et al. 2015;Zhou and Chen 2019).
In this study, XFEM based on Abaqus FEA software was used for the simulation of soil tensile behavior in order to investigate the role of crown cracks. The authors set a series of numerical simulations to explore the crown cracks formation and the effect of a tensile crack on soil slope stability. Factors, along with the position, strength and depth of cracks, were considered for sensitivity analysis. Combined with stress analysis, some interesting phenomena were found with conclusions made.

XFEM module
The Extended Finite Element Method (XFEM) is an improvement on FEM for the research of discontinuous processes, such as cracks. It is proposed by the Ted Belytschko's team (Moës et al. 1999). In the XFEM, a specific enrichment function is used for the discontinuity. When the enrichment function is applied to crack analysis, it can fit the asymptotic function of the tip and has a good expression for crack surface displacement jump. The standard XFEM approximation equivalent formula of u function in Abaqus is as follows: In Eq. (1), N I (x) is the node's normal shape function. The first term on the right of the equal sign is available for all nodes and is related to the continuous part of finite element. Displacement vector of the normal node is represented with u I . The second term applies to particular nodes, such as a shape function that supports a node cut inside a crack. Vector a I is the product of nodal enriched degree of freedom and H(x) denotes to the jump function across the discontinuous interface. The most restricted is third term, only for the shape function supporting the node cut off by the crack tip. The term is the product of nodal enriched degree of freedom vector b I . The elastic asymptotic properties of the crack tip are described with F (x). The jump function H(x) is described as follows: where, x represents the sampling point; x * is the point with the shortest distance from x on the crack and n stands for the unit vector of the crack outward normal at x * . Figure 2a shows the asymptotic function of crack tip in isotropic elastic material and the formula is as follows: In Eq. (3), ( r, ) is polar coordinate representation, whose physical meaning is that the origin is at the crack tip.
The node subset I * is the set of all nodes of elements cut by discontinuities. The global enrichment function can only function in elements whose nodes are all in the subset, I * . The level-set function is a scalar field whose zero-level represents discontinuity. The level-set function (x) , described as follows, determines whether an element is cut by discontinuities: In the above equation, I el represents the set of element nodes. The domain Ω is divided by the discontinuity into Ω P + and Ω P − and the level-set function can be positive or negative on either side of the discontinuity. In the domain Ω, the phantom node is used to describe the crack behavior and the node is initially superimposed on the real node previous to element separation. When there is no crack in elements, the phantom node corresponds to real nodes with complete constraint. If the element is divided into two parts by the discontinuity, corresponding phantom-nodes and real nodes will separate and no longer be tied together (Fig. 2b). They become interpolated by standard finite element shape functions as follows: where, h is the number of interpolated elements.

Tensile strength of soil
Tensile cracks are generated by soil stretching (Fig. 3a). A typical stress-displacement curve of the tension process is represented in Fig. 3b. Previous studies on tensile strength for soils (Hadas and Lennard, 1988) show formation of a tensile crack when the tensile stress ( t ) reaches the tensile strength ( f t ). Before t reaches f t , the t constantly increases with tensile displacement ( ▵ l ) (AB segment in Fig. 3b). In this stage, the tensile displacement consists of elasto-plastic deformation of soil ( ▵ l 1 + ▵ l 3 in Fig. 3a). When t reaches f t (point B in Fig. 3b), the soil is damaged and a crack is formed. This leads to stress accumulation at the disrupted part and the accumulated stress is continuously released along the crack, causing crack propagation and opening. In this stage, t decreases with the increase in ▵ l (BC segment in Fig. 3b), and the tensile deformation mainly consists of a crack opening, until it fully opens ( ▵ l 2 in Fig. 3a). At this time, the body of soil is completely separated (point C in Fig. 3b). Here, the authors define the length of the DC segment as failure displacement in the manuscript.
Abaqus has several traction separation laws for material damage and the maximum principal stress (MAXPS) failure criterion can be applied to the crack evolution of soil according to the analysis. In the MAXPS criterion, no crack is generated until the MAXPS reaches a certain value. Therefore, the value can be set as t for the soil. The criterion of MAXPS can be expressed as follows: where, o max is the maximum allowed principal stress, determined from the tensile strength of soil. Damage is initiated after the MAXPS reaches o max . The crack expands in the direction perpendicular to the MAXPS after the initial damage, and the evolution of an existing crack depends on the softening stage of the soil (BC segment in Fig. 3b). There are two types of ways to define the softening stage: defining . 3 Mechanical characteristics of the tensile crack of soil: a A sketch of tensile failure of materials, b an example of tensile stress-displacement curve of soil (modified from Tamrakar et al. 2005) displacement at failure or fracture energy. Displacement at failure is the length of DC segment, whereas the fracture energy can be considered as the area of DBC per length of material. The shape of softening phase curve (BC segment) can be determined by the form of index or discrete point data.
The composition of a soil slope usually includes clay, silt, sand and gravel, which are considered to be soil aggregates. To determine tensile parameters of soil aggregates, the authors studied some literature about tensile strength (Causarano et al. 1993;Hadas and Lennard 1988;Munkholm et al. 2002;Tamrakar et al. 2005;Zhang et al. 2006) and found that their tensile strength usually ranges from several kPa to tens of kPa for soil aggregates (Table 1). Many factors, such as density, water content, composition and porosity, affect the tensile properties of soils. The authors of this paper used the stress-strain curve, reported in Tamrakar's research (Tamrakar et al. 2005), as a typical stress-strain curve (Fig. 3b) and applied it to all the simulation where damage evolution was attenuated in a quadratic form.

Model validation
To evaluate whether XFEM can effectively analyze the tensile crack behavior of soil, a numerical tensile test was performed to calibrate it. The tensile test is referred from Tamrakar's study (Tamrakar et al. 2005). In this test, a compact clay-sand mixture specimen was stretched using a steel tensile mold at a steady velocity (Fig. 4a). The total mold surface area is 38.51 cm 2 with the total volume of 192.53 cm 3 . The minimum mold width is 3 cm and the depth is 5 cm. The soil specimen is prepared within this mold. Mechanical parameters of the steel mold taken as γ (density) = 7850 kg/m 3 , u (Poisson's ratio) = 0.29, and E (elasticity modulus) = 210 GPa. The soil specimen is Kanto loam, whose water content is 40% and dry density (ρ d ) is 0.7 g/cm 3 . The mechanical parameters of material were obtained from the research, whereas the damage parameters of soil were obtained from the stress-displacement curve ( o max = 10 kPa, the failure displacement = 0.2 mm). After the beginning of simulation, one end of the mold is fixed while the other moves at 0.35 mm/ min in the horizontal direction to pull the soil specimen, until the specimen is fractured and pulled into two parts. The results demonstrate the progress in crack propagation generated by stretching and the final crack figure in the simulation is consistent with the laboratory test ( Fig. 4b−e), indicating that XFEM can successfully simulate the tensile cracks in soil.

Slope stability analysis
In this section, a two-dimensional (2D) slope model was used to estimate whether XFEM can evaluate the slope stability with cracks. The model is an exemplary slope model with a toe angle of 45°, which has been widely used for the validation of slope stability (Bhandari et al. 2016;Wu et al. 2017). Its specific dimensions and boundary conditions are shown in Fig. 5a. Mechanical soil parameters were taken as γ = 2000 kg/m 3 , φ (internal angle of friction) = 34°, c (cohesion) = 10 kPa, u = 0.27, and E = 40 MPa from Bao's literature of soil aggregates (Bao et al. 2019). The tensile stress-strain curve was adopted from Tamrakar's research and corresponding soil damage parameters were taken in a relatively small value as o max = 1 kPa and the failure displacement = 0.02 mm. The stability of slope without XFEM element and the slope using XFEM with an internal crack were calculated.  Cai et al. (2017) In the simulation, the failure of slope relied on the strength reduction technique (Matsui and San 1992) and the strength reduction factor (SRF), which have been widely applied to determine the potential slip surface in slope (Niant et al. 2012;Jiang et al. 2015;Bao et al. 2020). The Shear Strength Reduction (SSR) technique discounts C and into C R and R using Stress Reduction Factor (SRF) and substitutes C R and R into the calculation until slope failure. At this moment, the discounted SRF is equal to Factor of Safety (FOS) of the slope. The SRR technique can be expressed as Eqs. (9) and (10) and it is performed by defining the field variables in this study: where, c is the cohesive force; φ is the internal friction angle; c R is the cohesive force after strength reduction; R is the internal friction angle after strength reduction; σ is the normal compressive stress on slip surface; τ is the shear strength per unit area on slip surface; and K is the strength reduction factor.  Fig. 4 Comparison result of real and simulation tensile tests: a-d Propagation of crack during the stretching based on the stress result, and e tensile test in Tamrakar's study (Tamrakar et al. 2005) The phenomenon that plastic zone extends from the toe to the upper of the slope was treated as the sign of failure (Shen and Karakus 2014). The Mohr-Coulomb (MC) failure criterion was adopted for the material. it can be considered that the cracks on slope might affect the shape of potential slip surface and even the value of FOS. In the following sections, we discuss several factors, including position, strength and depth of cracks that might influence the stability of slope based on the above model. Figure 6 shows the strength reduction of a slope without a pre-existing crack element when o max is equal to 5 kPa and the failure displacement is 0.1 mm. The simulation result is consistent with the field survey of most soil landslides. After strength reduction, the slope is deformed towards the free surface with slope toe swelling (Fig. 6a). Then, a part of the soil at the upper slope gradually changes from compressive stress state to tensile stress state (Fig. 6d). When the tensile MAXPS in the tension zone accumulates, damage starts ( Fig. 6e−f). The crack is formed at the slope crown and the soils presented on both sides of the crack have an obvious vertical displacement difference when the slip surface is completely penetrated (Fig. 6g−h). This presents a natural failure mode of soil slopes and proves that the crown crack is a signal of metastability.

Sensitivity analysis
It should be noted that the tensile crack appears and completely opens at the slope crown prior to slip surface penetration (Fig. 6e−f). This means that the length of a slip surface shearing is shorter than its entire length, because the crack occupies a certain range of it. If a model does not consider tension effect or tensile cracks (e.g., the conventional LEM assumes that the entire slip surface is generated by shear effect), the calculated FOS will be slightly larger than the real value. In the FEM simulation in our case, the FOS of the slope is 1.298 when the plastic zone is completely penetrated without the XFEM element, while the FOS is 1.295 when the plastic zone extends to the bottom of the tension crack with the XFEM element.
In addition, the authors found that the soil damage parameters have an influence on crack formation (Fig. 7). When the value of o max is less than 4 kPa, the crack forms inside the potential failure mass (Fig. 7a, c) and displacement of soil on both sides of crack is continuous at the time of slope failure (Fig. 7b, d). When the value of o max is more than 7.5 kPa, no crack is formed at the time of slope failure (Fig. 7g, h). Both situations are not consistent with field investigation, and seldom occur in reality (Tang et al. 2019). Thus, they are considered unreasonable. It is only then that the value of o max ranges from 4 kPa to 7.5 kPa, the simulation results are consistent with the landslide field study, as shown in Fig. 7e and Fig. 7f. This phenomenon can be attributed to soil strength. Bonds and friction exist between the soil particles. They are expressed as internal friction and cohesion in shearing, whereas they are expressed as tensile strength in tension. Therefore, the tensile strength of soil is not entirely independent of other strength parameters, such as shear strength. In the simulation, the value of cohesion is taken as 10 kPa, while the internal friction is taken as 34°, indicating that o max is a moderate value which is not too small or large. In other words, the tensile strength of soil in the scale of 4 kPa to 7.5 kPa is also moderate and appropriate. Although the relationship between shear and tensile strength of soil requires more study, this theory is supported by the results of numerical simulation. It provides a new calibration concept to roughly estimate the tensile strength of soil, especially for the field test which is demanding to conduct in a laboratory.
Besides the crown cracks, other types of tensile cracks, generated by cycles of wetting and drying (Konrad and Ayad 1997), weathering (Hales and Roering 2007) and desiccation (Peron et al. 2009) usually appear in a slope. To inspect whether tensile cracks affect 1 3 stability of slope, six groups of numerical tests were set for various situations, including the crack located at the position of trailing edge, inside the potential failure zone and out of the potential failure zone (Fig. 8). Considering that the strength of damage parameters might determine the outcome, five groups of strength data (the very low, low, medium, high and very high tensile strength) were used in the computation. The value of o max was taken as 0.5, 2, 5, 20 and 50 kPa and the corresponding failure displacement was taken as 0.01, 0.04, 0.1, 0.4 and 1 mm, according to Fig. 3b. The contour map of stress magnitude corresponding to Fig. 8a-f is shown in Fig. 9a- Figure 8a represents a crack at the slope crown. The result shows that it hardly has any impact on stability of slope, regardless of the length of crack and its strength. This is probably because the length of soil that can be sheared is very short at the slope crown or the pre-existing penetrating crack releases the accumulated stress and strain to the ground surface. Figure 8b shows the situation of a slope with a crack outside the potential failure zone. The results indicate that it is difficult for the tensile crack to influence the stability of slope, regardless of the length of crack and tensile strength of soil. The situation of a slope with a crack inside the potential failure zone is presented in Fig. 8c. The results show that, when the crack is very short and leaves some distance to the potential plastic zone, it does not affect slope failure, regardless of the strength of crack. When the crack is short, but located at the position of the potential plastic zone (Fig. 8d), it still does not affect the slope stability, regardless of the tensile strength of soil. This is due to shearing as the main effect at this part of the slip surface. Figure 8e shows a tensile crack, with relatively small damage parameters ( o max = 2 kPa), that is at the upper part of potential slip surface. Compared to Fig. 8d, the crack can be stretched to propagate to the ground surface after strength reduction. This is almost certainly because one tip of the crack is located at the tension zone (Fig. 6h), changing the shape of original slip surface and FOS in a smaller scale. Lower crack tip transmits stress, whereas the upper tip accumulates stress during strength reduction and it directs stress propagation. When the crack propagates to the ground surface, the accumulated stress is fully released (Fig. 9e). When the damage parameter ( o max = 50 kPa) is larger than the tensile MAXPS (13 kPa), the crack is not stretched, and the potential slip surface does not change. Figure 8f shows when a crack is much longer than the depth of original slip surface, the presence of a crack causes stress redistribution after strength Most importantly, the existence of cracks will cause discontinuities in stress (Fig. 9). A crack buried underground will transmit stress on one tip and accumulate stress on the other tip. If one crack tip propagates to the ground surface, all the accumulated stress will be released. Whether the tensile crack is pulled apart or not depends on the maximum tensile stress. According to simulation, with the exception of the tension zone on the upper slope, the slope is generally at a compressive state and it is hard for a short-medium tensile crack, in the compressive zone, to change the original stability. Only a long tensile crack, passed through the plastic zone, indicating a large scale of stress redistribution, will change the original stability.
To determine the stress variation in slope, from the initial deformation to failure, six nodes at different parts of the slope were set to monitor the corresponding information (Fig. 10a). The MAXPS, normal stress in X-direction (S11), normal stress in Y-direction (S22) and shear stress in XY plane (S12) were recorded. In Fig. 10b, SRF is less than 1.25, the MAXPS is less than zero and with the increase in SRF value the variation is very small. In this scale of SRF, the slope is at a compressive state and no tensile crack  is generated. When SRF is larger than 1.25, the MAXPS at the top of slope increases to a value of ~ 10 kPa, generating tensile cracks at the slope crown. Figure 10c shows the variation in S11 during the entire process. By comparing with Fig. 10b, it's observed that the direction of the MAXPS is generally horizontal, leading to vertical cracks. Figure 10d demonstrates that the slope of S22 is always less than a zero, due to the presence of gravity, and correlates positively with buried depth. In Fig. 10e, the shear stress continuously increases around the zone of slip surface, because the slip surface is mainly formed by the shear damage of soil. The stress state in slope illustrates the formation of crown cracks and why most cracks in Fig. 8 cannot determine the slope stability.

Conclusions
This study aims to illustrate the formation of crown cracks in a cohesive soil slope and to evaluate the effect of a tensile crack on slope stability using XFEM. The work is founded on a numerical slope model with a toe angle of 45° and cohesion of 10 kPa.  Fig. 10 Stress variation at different positions of slope: a Monitoring location of the model, b MAXPS of the monitoring nodes, c S11 of the monitoring nodes, d S22 of the monitoring nodes, and e S12 of the monitoring nodes with different damage parameters or lengths are set at different positions in a slope model to achieve the goal. Some conclusions are drawn from the simulation: 1. The MAXPS based on XFEM efficiently simulates the tensile behavior of soil. The cracks will be located according to the stress field and can develop within the computation. The FOS of the slope, considering tensile cracks, is slightly smaller than the FOS of the slope without considering tensile cracks; 2. With slope deformation, the tension zone gradually appears on top of slope. When crown cracks form, the potential slip surface is nearly to be penetrated; and. 3. Most of the tensile cracks in slopes hardly influence the slope stability, excluding the cracks located in the tension zone or cracks passing through the potential plastic zone with a long length.
In some aspects, these conclusions can guide the practical engineering. Namely, it is challenging to determine the tensile strength of soil slopes. However, researchers can obtain the soil tensile parameters by calibrating the position of crown cracks, based on the numerical simulation of slope stability assessment. They can also estimate the evolution of slope deformation by comparing the crown cracks depth or the depth of drop heads between the actual slope and numerical model.
Previous conclusions are based on the condition that seepage in cracks is not considered. In fact, water easily penetrates the crack , causing propagation (Zhou et al. 2009a, b) and even generating landslide (Zhang et al. 2012;Xing et al. 2021). In addition, the water flow along the crack will cause an increase in pore water pressure in the lower soil and produce lubrication, macroscopically demonstrating the reduced shear strength in lower soil. Nevertheless, the ways and means to consider these complex interactions and behaviors in a numerical model is still a challenging task.