Formation mechanism and stability analysis of the DR landslide in national road 534 of Datian, Sanming, Fujian Province, China

Abstract Rainfall-induced landslides occur frequently along roads in Fujian Province, China. Previously, landslides located at roads in Jianou and Yanping in Fujian were reported to claim 23 and 24 lives, respectively. The DR landslide (located on National Road 534 of DR village, Datian County, Sanming City, Fujian Province), with a state of instability and the potential of burying the road and buildings and causing human casualties, is representative. Therefore, there is a need to analyze the stability of the DR landslide. Through detailed reconnaissance and modelling, this study investigated the characteristics, formation mechanism and stability of the DR landslide. Based on GeoStudio software, the stabilities of the upper, middle and lower parts of the landslide were explored in natural, rainfall, earthquake and coupled rainfall-earthquake conditions, with the finite element, Bishop’s simplified methods and others. The results demonstrate that under the above four conditions, the upper and middle parts were stable with stability coefficients greater than the standard engineering values for China, while the lower part was unstable with stability coefficients of critical slip surfaces of 0.871, 0.798, 0.867 and 0.727, respectively, which are less than the standard values. Compared with the natural condition, rainfall influenced stability less, whereas the earthquake and coupled conditions greatly affected the stability. Based on the analyses, some technical support is provided to implement landslide prevention, mitigation and control measures such as filling cracks with grouting and widening the retaining wall.

1. Introduction method, in particular, is an acknowledged method that provides sufficient accuracy and meets engineering requirements by considering the horizontal forces between soil slices while ignoring the friction between them (Writing Group of Professional Voluntary Standard of the People's Republic of China 2019). The elements used to assess slope stability under earthquake action consist mainly of permanent displacement and stability coefficient. The accelerations at which the block average accelerations exceed the yield acceleration are employed to obtain the permanent displacement by quadratic integration with the precondition of a landslide body as a rigid-plastic sliding mass, that is, the Newmark method (Newmark 1965). Furthermore, Liu et al. (2003) presented the minimum mean stability factor (the calculation equation can be found in Section 2.2.2) to evaluate slope stability under the action of an earthquake because the minimum stability coefficient arose at one moment during the process of an earthquake and was not appropriate for analyzing stability. Currently, theory calculation methods have been gradually developed and applied to evaluate slope stability (Chi et al. 2006;Qin et al. 2006;Lin 2008;Zhao and Liao 2012;Liu et al. 2018;Wang et al. 2018;Vadivel and Sennimalai 2019). However, by using theory calculation methods, the megascopic deformation and distribution of stress and strain in slopes are invisible and cannot be explored.
Physical model tests can overcome the drawbacks of theory calculation methods. Montrasio et al. (2016) performed a flume test to simulate the triggering process of shallow landslides induced by rainfall. Zhang et al. (2017) investigated the seismic responses of a loess-mudstone landslide model at various earthquake wave amplitudes based on centrifuge shaking table tests. Shen et al. (2020) studied the evolution process of the Dehua Shishan landslide located in Fujian Province under rainstorm action using a physical model test. However, such tests tend to require extensive resources, including time, effort and money. Compared with physical model tests, numerical simulations not only produce results simply and quickly by computer calculation, but can simulate more complex processes of landslides under various conditions, which have been widely used in the analysis of various aspects of landslides (Collins and Znidarcic 2004;Kim et al. 2004;Qi and Vanapalli 2015;Wolter et al. 2016;Chen et al. 2020). Common numerical simulation methods mainly include finite element, finite difference and distinct element methods. For example, under rainfall and earthquake, Cui et al. (2008) analyzed the stability of Badashan landslide group located in Qinghai Province in SEEP/W and SLOPE/W modules of GeoStudio software and its deformation characteristics in FLAC3D finite difference software; Gischig et al. (2015) explored interactions of deep-seated rock slope instabilities and seismic waves by mainly using UEDC distinct element software; Acharya et al. (2016) carried out seepage and stability simulations of seven slopes under rainfall based on SEEP/W and SLOPE/W modules of GeoStudio software. However, the results of numerical simulations may cause significant differences compared with practical situations. Therefore, model verification is essential for simulation results in accordance with practical situations. To address the aforementioned shortcomings, displacement back analysis of engineering geology, that is, observed displacements in the numerical model correspond to site displacements by tuning material parameters, should be adopted. Several landslide precautions have been mentioned in previous literature, such as retaining wall support, bolt anchorages and plant reinforcement technology (Cui et al. 2008;Lv et al. 2017;Chen et al. 2020). In summary, the aforementioned methods and controlling measures provide technical support and valuable references for studying the DR landslide.
Based on the literature, theory calculation methods and numerical simulations were used in the present study. These not only calculate slope stability factors, but also display megascopic deformations and stress and strain distributions of slopes, and are feasible, sophisticated and advanced methods. Furthermore, rainfall and earthquakes are crucial factors that affect slope stability (Collins and Znidarcic 2004;Kim et al. 2004;Cui et al. 2008;Gischig et al. 2015;Qi and Vanapalli 2015;Acharya et al. 2016;Wolter et al. 2016;Lv et al. 2017;Wu et al. 2017;Chen et al. 2020), which indicates that exploring the influence of rainfall and earthquakes on slope stability is highly significant. A typical profile (1-1 profile) of the DR landslide, divided into three regions, that is, the regions I, II and III (the upper, middle and lower parts, respectively), was utilized to analyze the stability and characteristics of deformation and failure of the three regions under four simulated conditions (natural, rainfall, earthquake and coupled rainfall-earthquake conditions). Based on the SLOPE/W, SEEP/W, SIGMA/W and QUAKE/W modules in GeoStudio software, these were performed using Bishop's simplified, finite element, Newmark and minimum mean stability factor methods. Finally, a displacement back analysis of engineering geology was employed for model verification.  and topography increase from high in the northwest to low in the southeast. Most elevations range from 370 to 710 m a.s.l., with the highest reaching 1369 m a.s.l. There is a complex terrain in the slope front and a valley-surrounding-shaped depression with steep terrain on the back wall.

Data and methods
2.1.1.2 Lithology. According to the drilled samples, the following strata are chronologically listed: (1) the plain fill of Quaternary Holocene (Q 4 ml ) consists of loose cohesive soil generated by the landslide, mixed with 10% sand and a small amount of gravel; (2) silty clay and silty clay containing gravel and rubble constitute Quaternary Holocene (Q 4 dl ) slope deposits that mainly include silt and clay. The former and latter contain quartz sand and 20-35% gravel, respectively; (3) as the subterrane, the Tianwadong Formation of the Upper Devonian series (D 3 t) comprises completely weathered glutenite, sand-soil-shaped strongly weathered glutenite, sand-soil-shaped strongly weathered sandy mudstone, block-shaped strongly weathered glutenite, moderately weathered glutenite and slightly weathered glutenite. The attitude of the rocks is 200 /29 . 2.1.1.3 Tectonic structure. The tectonic unit of this area is located in the depression belt of southwest Fujian. According to the Seismic Ground Motion Parameters Zonation Map of China (GB 18306-2015) (China Earthquake Administration 2016), the basic seismic intensity is 6 degrees; the horizontal peak ground acceleration (PGA) is 0.05 g; the characteristic period of response spectra is 0.35 s.
2.1.1.4 Hydrogeological conditions. The developed groundwater, mainly recharged by atmospheric precipitation, is comprised of phreatic pore water in Quaternary slope sediment and pore-fissure water in weathered bedrock. Water volume is significantly affected by rainfall. Depending on the location, the depth of groundwater was 1.00 to 19.10 m. The average depths of groundwater were 3.70, 5.30 and 7.90 m in the regions I, II and III, respectively. In terms of the 1-1 profile, the groundwater depths were 4.7, 5.0, 1.00, 3.00 and 5.6 m in drill holes Nos. 1, 2, 3, 4, and 5, respectively ( Figure 2).   c and c' are natural and saturated cohesion, respectively; U and U' are natural and saturated internal friction angles, respectively; n is the porosity; provide the following specifications for evaluating the stability of a landslide. Under natural conditions, when the stability factor (F s ) is less than 1.0, the landslide is in an unstable state; when F s is greater than or equal to 1.0 and less than 1.05, the landslide is in a metastable state; when F s is greater than or equal to 1.05 and less than the factor of safety, the landslide is basically stable; when F s is greater than or equal to the factor of safety, the landslide is stable. In addition, the safety factor has a large value of 1.30, in that the landslide has caused serious harm to the road. Under rainfall or earthquake conditions, the stability coefficient of the landslide should not be less than  year. May and June constitute the primary flood season in Sanming. The rock-soil mass gradually became saturated, and the sliding force increased due to the influence of persistent heavy rain during May and June 2019, causing further sliding. The partially broken anchor heads of the anti-sliding piles and permeable walls were ascribed to bulging cracks in the protective face and retaining wall at the toe and tensile cracks in the rear ( Figure 3). As a result, the landslide has posed a serious threat to human life, buildings, and National Road 534. Based on detailed reconnaissance and other data, the main sliding direction is approximately 170 with an inclination from 15 to 20 . The landslide has turned into a cemetery with nearby lush vegetation. A protective face and retaining wall lie at the toe of the slope, where a meandering road exists. An overview of the landslide is shown in Figure 4. It was unstable with a buried depth of overall slip surfaces of 3-22 m. The total length and width were approximately 236 m and 58-137 m, respectively. In the region I, the buried depth of the slip surfaces was 3-8 m, the total length and width of the sliding mass were approximately 66 and 58 m, respectively, and the volume of the sliding mass was 1.914 Â 10 4 m 3 . In the region II, the buried depth of the slip surface was 3-12 m, the total length and width of the sliding mass were approximately 63 and 65 m, respectively, and the volume of the sliding mass was 4.095 Â 10 4 m 3 . In the region III, the buried depth of the slip surfaces was 6-22 m, the total length and width of the sliding mass were approximately 107 and 137 m, respectively, and the volume of the sliding mass was 17.591 Â 10 4 m 3 . The 1-1 profile ( Figure 2) shows that the elevations of the top left, top right, and bottom are approximately 595, 517, and 480 m a.s.l., respectively. The length of the bottom is 300 m. There are five slip surfaces numbered Nos. 1, 2, 3, 4, and 5, from the region I to the region III, in which Nos. 1 and 2 are located in the region I, No. 3 lies in the region II, and Nos. 4 and 5 are situated in the region III.
According to interviews with local villagers, the personnel in one reconnaissance institute of Fujian analyzed the slope stability by theoretical calculation.

Numerical route
The 1-1 profile created in the AutoCAD software was imported into the GeoStudio software for modelling. In the model, the heights of the left and right boundaries are 114.58 and 37.24 m, respectively. The length of the bottom is 300 m ( Figure 5). First, the existing circular slip surfaces were specified to solve in GeoStudio software. Next, the range of the entry and exit of slip surfaces was specified so that slip surfaces were searched automatically and then solved again. Finally, the critical slip surface was identified in the most dangerous region. On the basis of the above, the stability of the landslide was investigated under natural (original state), rainfall (current state), earthquake (future state), and coupled rainfall-earthquake conditions (future state). The original slope was reproduced under a natural condition. The slope during the current state was modelled under a rainfall condition. The slope during the future state was simulated in earthquake and coupled rainfall-earthquake conditions. Under the natural condition, the stability of the slope with the natural water table was analyzed using the SLOPE/W module. Under the rainfall condition, the steady-state and transient seepage results were first produced in the SEEP/W module, and then the results were imported into the SIGMA/W and SLOPE/W modules for slope deformation and stability analysis, respectively. In earthquake conditions, the initial static and dynamic analyses were first performed in the QUAKE/W module, and then the results served as the initial conditions in the SLOPE/W module for the stability calculation. In coupled rainfall-earthquake conditions, the calculated results in the SEEP/ W module acted as the initial condition in the QUAKE/W module, and the obtained results were then introduced into the SLOPE/W module to investigate the slope stability.
Bishop's simplified method, based on the limit equilibrium theory, was employed to calculate the stability factors of the landslide under natural and rainfall conditions. It is expressed as Eq. (1) (Qian 1988): where W i is the weight of slice i, b i is the width of slice i, l i is the pore water pressure of a concentrated point of the slice base, h i is the slice base inclination, and c' and u' are the indices of the effective stress strength of the landslide mass. Under earthquake action, the minimum mean stability factor method and Newmark method were adopted. Liu et al. (2003) presented the minimum mean stability factor, expressed as Eq. (2), to evaluate slope stability in the event of an earthquake: where F smin is the minimum mean stability coefficient, F s0 is the stability coefficient in a static state, and F smin is the minimum stability coefficient with variation in ground motion. The Newmark method is unable to determine whether a slope is stable or unstable because of the absence of a standard of damage (Zheng et al. 2010). Therefore, by combining the Newmark method with the minimum mean stability factor method, the stability of the slope was finally assessed using the minimum mean stability factor as an index in both earthquake and coupled rainfall-earthquake conditions.

Formation mechanism
Factors affecting the stability of the DR landslide mainly consist of human engineering activities, rainfall, an adverse terrain environment, lithology and the structure of the rock-soil mass, rendering the landslide susceptible to failure. The path of Provincial Road 306 was altered, resulting in the existence of 20-25 m of sloped surface and a free face, which disrupted the initial stress balance within the slope. Furthermore, the U-shaped topography of the slope favours the accumulation of groundwater in the sliding mass. In the region III, the layer was comprised of three thick layers, one soil layer, another loose and completely weathered layer and the other strongly weathered layer; they were easily softened by water. In addition, the moderately weathered layer below served as a relative aquifuge. As a result, the front of the slope body was readily saturated. Similarly, in the regions I and II, the rocksoil mass became easily saturated because of the readily softened upper soil and moderately weathered layer below serving as a relative aquifuge. The bedding plane and slope surface were oriented in the same direction, with an inclination of 15 to 20 . Heavy rain infiltration, saturated rock-soil mass, increased body weight and the decreased shear strength all contributed to the instability of the slope body. In summary, the landslide is unstable because of the combined effects of human engineering activities, rainfall, an adverse terrain environment, lithology, and the structure of rock-soil mass.

Natural condition
Under natural conditions, the SLOPE/W module in GeoStudio software was employed for the stability analysis. The Mohr-Coulomb plasticity model served as the constitutive model of the rock-soil body. A water table was then drawn. The model is illustrated in Figure 5. The calculated results of Bishop's simplified method are shown in Figure 5. In the region I, the stability factors of slip surfaces Nos. 1 and 2 were both greater than 1.30, indicating that the region I was stable. In the region II, the stability coefficient obtained from slip surface No. 3 was greater than 1.30, demonstrating that the region II was stable. In the region III, the stability factors obtained from slip surfaces Nos. 4 and 5 were less than 1.30. However, the stability coefficient of the critical slip surface located in the region III was less than 1.0 (Figure 5b), indicating that the region III was unstable. That is, the landslide was most likely to occur in the region III.

Rainfall condition
The reconnaissance report indicated that the variation in the water table in the slope was approximately 2 m. Rainfall conditions were defined as the water table rising by 2 m cumulatively in eight days after considering the rainfall history in Sanming City.
The changes in steady-state and transient seepage were examined using the SEEP/W module in the GeoStudio software. To investigate transient seepage, eight time steps were adopted. The mesh was comprised of quads and triangles of a finite element mesh divided into 1595 nodes and 1828 elements. The rock-soil masses were unsaturated and saturated above and below the water level, respectively. The soil water characteristic curves and hydraulic conductivity functions of plain fill, silty clay, and silty clay containing gravel and rubble were acquired using sample functions (water content functions of different types of soils) and the van Genuchten method (van Genuchten 1980) in GeoStudio software, respectively (Figures 6 and 7). The van Genuchten method is expressed as Eq. (3): where k w is the soil hydraulic conductivity, k s is the saturated hydraulic conductivity, a, m, n are curve fitting parameters, among them, n ¼ 1=ð1 À mÞ, and w is the required matrix suction range. As for boundary conditions, the zero-flux boundary was utilized on both sides above the water level; the fixed head boundary was used on both sides below the water level, and the zero pressure head boundary was adopted on the water level.
As the water table rises, the unsaturated rock-soil mass gradually becomes saturated, and the pore water pressure increases. When the pore water pressure near slip surfaces increased, the effective stress and shear strength of slip surfaces tended to decrease, which intensified the impending slope failure. Therefore, it was necessary to explore the effects of pore water pressure on the landslide stability. According to the SEEP/W module solution, whether in steady-state or transient seepage cases, the maximum positive pore water pressure emerges in the deeper parts of the slope (below slightly weathered rock), while the minimum negative pore water pressure arises on the slope surface. The seepage field varied slightly. In the initial state, the maximum and minimum pore water pressures were 1069.8 and -92.271 kPa, respectively (Figure 8a). After the water level rose to 2 m cumulatively in 8 days, the maximum pore water pressure was 1089.4 kPa, which was 1.83% higher than that in the initial state, whereas the minimum pore water pressure remained unchanged (Figure 8b). Some points were selected as sites for monitoring the pore water pressure near the slip surfaces (Figure 9). Note that the pore water pressure near the slip surfaces increased over time, but the overall growth was minor ( Figure 10). With respect to the aforementioned analyses, the growth of pore water pressure had a relatively less influence on inducing further sliding. The pore water pressure produced in the SEEP/W module was taken as the initial pore water pressure in the SIGMA/W module to analyze slope deformation. In the SIGMA/W module, the mesh was the same as that in the SEEP/W module. A linear elastic constitutive model was used in this study. The boundary conditions along the model laterals were x-direction freed and y-direction fixed; the boundary condition along the model bottom was the xand y-directions fixed. Figure 11 shows the distribution of the horizontal and vertical displacements. The horizontal and vertical displacements ranged from -0.15 to 15.61 cm and -31.4 to 0 cm, respectively. Both peaks appeared at the outer edge of region III. The outer edge of region III was the main gathering location of the horizontal and vertical displacements, where there was a concentration of greater displacements. In the external and rear edges of region I and the outer edge of region II, both displacements were minor. Consequently,  deformation appeared in the slope under this circumstance, especially serious deformation in the region III.
The final results were produced in the SLOPE/W module by importing the calculated results in the SEEP/W module (Table 2 and Figure 12). On the first day, the stability coefficients obtained from slip surfaces Nos. 4 and 5 were slightly greater than those under the natural condition. Subsequently, they gradually declined. It was projected that the water had pressed against the front of the slope due to a period of water permeation, causing the anti-sliding force to increase temporarily. The slope stability gradually deteriorated with an increase in the water table. As a result, the stability coefficients obtained from slip surfaces Nos. 4 and 5 first increased and then  declined. On the eighth day, the stability factors of slip surfaces Nos. 1 and 2 in the region I were greater than 1.20, indicating that the region I was stable. The stability coefficient of slip surface No. 3 in the region II was greater than 1.20, indicating that the region II was stable. The stability coefficients obtained from slip surfaces No. 4 and 5 in the region III were less than 1.20 and the stability factor of the critical surface located in the region III was less than 1.20, demonstrating that the region III was unstable. Compared with those under the natural condition, the stability coefficients of the critical slip surface fell by 2.2%, 4.0%, 5.2%, 6.1%, 6.8%, 7.3%, 7.9% and 8.4% from the first day to the eighth day, respectively. With respect to the preceding analyses, the regions I and II exhibited smaller deformation, while the region III exhibited larger deformation. The increase of pore water pressure had relatively less influence on further sliding, but the defined water level promoted landslide instability, that is, the stability decreased over time. If rainfall continues for a long time, slope failure will be inevitable.

Earthquake condition
The basic seismic intensity of the area is six degrees. Therefore, the present study only considered horizontal seismic acceleration with a peak ground acceleration (PGA) of 0.05 g. A new seismic wave record, derived from correcting the peak ground acceleration of an existing typical seismic wave record, was adopted owing to the lack of local seismic wave records. The duration of the seismic wave was 50 s with a time interval of 0.05 s. The PGA was -0.05 g at 16.865 s, as illustrated at point A in Figure 13. The mesh in the QUAKE/W module was identical to that in the SEEP/W module. A linear elastic material model was used in this study. In the initial static analysis, the boundary conditions were x-direction fixed and y-direction freed along the model laterals and the xand y-directions fixed along the model bottom. In the dynamic analysis, the boundary conditions were x-direction freed and y-direction fixed along the model laterals and the xand y-directions fixed along the model bottom. Additionally, history points 1, 2, and 3 were located in the regions I, II, and III, respectively, and were used to record the displacements (Figure 14).
In the dynamic response analysis, the displacements of the three monitoring points were obtained. Because only the horizontal seismic force was considered, the variation in the horizontal displacement was examined; that is, the variation of the vertical displacement was ignored. The findings revealed that the variation trends of the displacements were approximately identical ( Figure 15). The displacement variation at history point 3 is as follows: based on the modelling, from the start to the 5.765th second, the positive growth of displacements proceeded slowly. From the 5.765th second to the 11.33th second, the displacements fluctuated greatly and grew markedly in the negative direction. From the 11.33rd second to the 21.73rd second, undergoing the action of the peak ground acceleration at 16.865 s, the displacement trend exhibited a positive growth rate with smaller rates in contrast with those from the previous period. From the 21.73rd second to 50th second, that is, the duration corresponding to the curve segment of the horizontal earthquake acceleration record (Figure 13), the displacements increased more slowly in the negative direction because of the minor earthquake intensity. The other two history points had similar displacement variations. Finally, the displacements of the three monitoring points were approximately -0.15 m.
After a simulated 50-second earthquake, the slope deformation was concentrated on the outer and rear edges of region I and the outer edge of region III, that is, the deformation was mainly focused on plain fill, silty clay and silty clay containing gravel and rubble, and less on the bedrock below ( Figure 16). Under the horizontal earthquake load, the overall horizontal peak displacement was approximately 14.7 cm, with the maximum displacement at the outer edge of region I (Figure 17a). Moreover, the outer edges of all regions were the main gathering locations of the   displacement, where the displacement increased from the inside to the outside and was more concentrated. In addition, the vertical peak displacement was 0-0.31 cm, with the maximum displacement occurring at the outer edge of region I (Figure 17b). Furthermore, the displacement was concentrated mostly on the external and rear edges of region I and the outer edge of region III, where it rose from the inside to  the outside and possessed a greater concentration. Concerning the horizontal and vertical peak displacements, greater concentrations appeared in the external edges of regions I and III, indicating that these regions were susceptible to failure.
The strain field analysis ( Figure 18) showed that the maximum shear strain was 4.22 Â 10 À9 through 8.29 Â 10 À5 . The maximum shear strain was distributed on the external edges of all regions, especially the outer and rear edges of region I. The maximum was observed in a small region near the slope crown, with less effect on the entire slope. Consequently, the modelling revealed that the region I was prone to failure.
The final results, obtained by combining the time-history curve of the stability coefficients with Eq. (2), are shown in Figure 19 and Table 3. The region I was stable because the minimum mean stability factors of slip surfaces Nos. 1 and 2 in the region I were greater than 1.15. The minimum mean stability coefficient of No. 3 slip surface in the region II was greater than 1.15, indicating that the region II was stable. The minimum mean stability coefficients of slip surfaces Nos. 4 and 5 and the critical slip surface located in the region III were less than 1.15. In other words, the region III was unstable.
According to the above results, the regions I and III are believed to be susceptible to failure based on the displacement field analysis. In addition, the region I was subject to failure based on the strain field analysis. Overall, the stability coefficients under the earthquake condition were quite different from those under the natural condition. Hence, the earthquake would have a greater impact on the stability of the DR landslide.

Coupled rainfall-earthquake conditions
The pore water pressure obtained from the SEEP/W module acted as the initial pore water pressure in the QUAKE/W module, and the calculated results were imported into the SLOPE/W module to study the stability of the landslide under the coupled conditions of rainfall and earthquake.
Under coupled conditions, the variation of pore water pressure was identical to that under the rainfall condition, and the variation trends of the horizontal Figure 19. Time-history curves of stability coefficients of slip surfaces in the earthquake condition.
displacements of the three monitoring points were very similar to those under the earthquake condition, with only minor differences ( Figure 20). However, the displacement and strain fields were different from those under the earthquake condition.
Under the coupled rainfall-earthquake actions, the deformation was focused on the outer edges of all regions and the rear edge of region I. In other words, the deformation was concentrated primarily on plain fill, silty clay and silty clay containing gravel and rubble and less on the bedrock below. In contrast to the earthquake condition, a more obvious deformation appeared in the region II ( Figure 21). Under the coupled conditions, the overall value of horizontal peak displacement was approximately 14.7 cm, with the external edge of region I exhibiting the greatest displacement ( Figure 22a). Moreover, the concentrated displacement was mainly located at the outer edges of all regions, where the displacement increased from the inside to the outside, and a greater concentration of displacement occurred. Additionally, the value of the vertical peak displacement was 0-0.24 cm (Figure 22b). Its peak was observed in the outer and rear edges of region I and the outer edge of region III. Moreover, the displacement was primarily focused on the external and rear edges of region I and the outer edge of region III, where it increased from the inside to the outside, and a greater concentration of displacement emerged. Overall, the horizontal displacement was greater and the vertical displacement was smaller in comparison with those under the rainfall condition. Presumably, the reason for this was that the shearing action of the horizontal earthquake wave may suppress the body downhill. Note: F s0 is the stability coefficient in a static state; F smin is the minimum stability coefficient with the variation of ground motion; F smin is the minimum mean stability coefficient. Figure 20. Time-history curves of horizontal displacement at history points in coupled rainfallearthquake conditions.
Regarding the horizontal and vertical peak displacements, a greater concentration of both displacements was observed in the outer edges of regions I and III, indicating that these regions were subject to failure. In the strain field analysis (Figure 23), the value of maximum shear strain ranged from 4.38 Â 10 À9 to 5.80 Â 10 À5 . The maximum shear strain was focused on the  external edges of all regions, especially on the outer and rear edges of region I and the external edge of region II. The peak appears in a small region near the slope crown and in the outer edge of region I, with less effect on the entire slope. As a result, the regions I and II were susceptible to failure.
The minimum mean stability factors were derived by combining the time-history curve of the stability coefficients with Eq. (2). The final results are shown in Figure 24 and Table 4). The region I was stable because the minimum mean stability factors of slip surfaces Nos. 1 and 2 were greater than 1.20. In the region II, the minimum mean   Note: F s0 is the stability coefficient in a static state; F smin is the minimum stability coefficient with the variation of ground motion; F smin is the minimum mean stability coefficient.
stability coefficient of slip surface No. 3 was greater than 1.20, indicating that the region II was stable. In the region III, the minimum mean stability factors of slip surfaces Nos. 4 and 5 were less than 1.20. Furthermore, the minimum mean stability coefficient of the critical slip surface located in the region III was also less than 1.20, indicating that the region III was unstable. Based on the aforementioned results, the regions I and III were vulnerable to landslide failure in the displacement field study; the regions I and II were prone to failure in the strain field analysis. Overall, the stability coefficients under the coupled conditions were more different than those under the rainfall condition, but they were less different than those under the earthquake condition. Hence, in terms of the stability of the landslide in coupled rainfall-earthquake conditions, the earthquake would be expected to play an important role, whereas rainfall would only play a minor role. Table 5 displays the results calculated for the stability of the DR landslide. In all modelled conditions, the analyses indicated that the regions I and II were stable, while the region III was unstable. The rainfall, earthquake, and coupled rainfall-earthquake actions all contributed to instability. The stability factors of the critical slip surface under rainfall, earthquake and coupled rainfall-earthquake conditions fell by 8.38%, 0.45% and 16.54%, respectively, compared with those under the natural condition. The results exhibited a monotonically declining trend in the stability factors for slip surfaces Nos. 1, 2 and 3 under the four conditions. This trend, however, did not appear in Nos. 4 and 5, and the critical sliding surfaces. Slip surfaces Nos. 4 and 5, with a relatively greater range, were found in a gentle terrain favouring water accumulation. In coupled conditions, we inferred that the press action of water and shearing action of the horizontal earthquake wave may make Nos. 4 and 5 sliding surfaces more stable compared with the earthquake condition. However, the critical sliding surfaces, with a relatively smaller range, lying in a steep terrain that does not favour water accumulation, may become more unstable. Furthermore, Eq. (2) is used to calculate the minimum mean stability factors, causing the error, probably due to the fluctuation of stability factors under earthquake.

Comprehensive analysis
Overall, in contrast with those under natural conditions, the stability coefficients under the rainfall condition had few differences, while the stability coefficients in Note: F s is the stability coefficient in natural and rainfall conditions; F smin is the minimum mean stability coefficient in earthquake and coupled rainfall-earthquake conditions. earthquake or coupled rainfall-earthquake conditions were significantly different. Consequently, the rainfall action had a small effect on the stability of the landslide, while earthquake or coupled rainfall-earthquake action had a greater impact on the landslide stability. In contrast, the landslide was more likely to fail under coupled action.

Discussion
Based on the results of the numerical modelling, the stability of the DR landslide under the natural condition and the influence of rainfall and earthquake on its stability and characteristics of its deformation and failure were discussed. The results showed that the regions I and II were stable, while the region III was unstable under natural (original state), rainfall (current state), earthquake (future state), and coupled rainfall-earthquake conditions (future state). Under the rainfall condition, the regions I and II were stable with smaller deformation, while the region III was unstable with larger deformation, which corresponded with the current situation that bulging cracks in the protective face and retaining wall at the toe and tensile cracks in the rear were observed. This indicates that the results are reliable. The DR landslide was unstable under the natural condition. By comparison, its stability was worse under rainfall, earthquake, and coupled rainfall-earthquake conditions, which conformed to the facts. This indicates that the results obtained are trustworthy.
According to interviews with local villagers, the personnel in one reconnaissance institute of Fujian used a theoretical calculation method to analyze the stability of the DR landslide. Less control of the landslide was conducted. The landslide is likely to slide further on a larger scale, posing severe hazards such as damage to the cemetery, buildings, the road at the toe of the slope, and human casualties. Therefore, there is an urgent need to gain insights into its stability. Previous studies on landslides along roads in Fujian have mostly focused on the use of theoretical calculations for preliminary analyses of landslide stability (Liao 2003;Chi et al. 2006;Huang et al. 2006;Qin et al. 2006;Lin 2008;Wu and Huang 2011;Ye 2012;Zhao and Liao 2012). However, the stability and features of deformation and failure of the DR landslide under various conditions were comprehensively investigated by combining theoretical calculation methods with the finite element method, which can enhance the study level of the area. Not only do the aforementioned findings have referential value so that they can be used to carry out research in the area, but they also offer technical support to implement prevention, mitigation and control measures for the landslide.
Some noteworthy points are as detailed below.
1. In terms of the landslide simulated with GeoStudio software, further deformation and failure along cracks were not taken into account, and the software did not allow consideration of sliding behaviour if the slip mass separated from the slip bed. Therefore, the slope failure process needs to be further explored. 2. The DR landslide is still developing at the time of writing. If the area is subject to rainfall, earthquakes, or other events, the slip surfaces may be activated, causing the region III to slide along slip surfaces; if this occurs, as expected, the slide masses in the regions I and II may also follow the slope downhill. A future landslide may result in severe hazards, such as damage to the cemetery, buildings, and the road at the toe of the slope, as well as human casualties. At present, prevention, mitigation, and control measures are needed to control the landslide, including filling cracks with grouting, widening the retaining wall, and adopting rock bolts supported at the toe or in the retaining wall. Further action should be taken to verify the effectiveness of the numerical simulation by monitoring displacements, groundwater, the physical field, earthquakes, and precipitation, among others.

Conclusions
The formation mechanism of the DR landslide was summarized according to the detailed reconnaissance. Based on the Geostudio software, the stability and characteristics of deformation and failure of the DR landslide were analyzed under natural, rainfall, earthquake, and coupled rainfall-earthquake conditions using Bishop's simplified, finite element, Newmark and minimum mean stability factor methods. In light of these analyses, the following points summarize our main conclusions: 1. The DR landslide was found to be unstable due to the interaction of human engineering activities, rainfall, an adverse terrain environment, lithology, the structure of rock and soil mass, and other factors. 2. Under the natural condition, the stability coefficients of regions I and II of the DR landslide are greater than 1.30, while the stability factor of the critical slip surface of region III is 0.871, which is less than 1.0. Under the rainfall condition, the stability coefficients of regions I and II are greater than 1.20, whereas the stability factor of the critical slip surface of region III is 0.798, which is less than 1.20. Under the earthquake condition, the minimum mean stability factors of regions I and II are greater than 1.15, while the minimum mean stability factor of the critical surface of region III is 0.867, which is less than 1.15. Furthermore, under coupled rainfall-earthquake conditions, the minimum mean stability factors of regions I and II are greater than 1.20, while the minimum mean stability coefficient of the critical surface of region III is 0.727, which is less than 1.20. Consequently, under the above conditions, the stability coefficients of regions I and II are greater than the standard values, i.e. the regions I and II are stable, while the stability factors of critical slip surfaces situated in the region III are less than the standard values, that is, the region III is unstable. Overall, in contrast to the stability coefficients under the natural condition, rainfall has a small effect on the stability of the landslide, while earthquake and coupled rainfall-earthquake actions have a significant influence on the stability. Moreover, the DR landslide is believed to be prone to failure under the coupled action described above. 3. Under rainfall conditions, the increased pore water pressure promotes slight sliding, however, the landslide stability decreases over time. The regions I and II with smaller deformation are in a stable state, while the region III with larger deformation is in an unstable state, which is in accordance with the current status of the DR landslide. Reinforcement measures and displacement monitoring must be implemented. 4. In earthquake conditions, deformation is predicted to occur at the outer and rear edges of region I and the external edge of region III. The deformation is projected to be concentrated on the outer edges of all regions and the rear edge of region I under coupled rainfall-earthquake conditions. According to the displacement field analyses, the regions I and III with a concentration of greater displacement will be vulnerable to failure under the above two conditions. Additionally, based on the strain field analyses, the region I with the occurrence of a concentration of greater strain will be prone to failure in the earthquake condition, and the regions I and II with a concentration of greater strain will be susceptible to failure in coupled rainfall-earthquake conditions.