Section 2

## Geomorphic Hazards

#### **Chapter 2**

### Genetic Algorithm Based Software for Optimization and Design of Piles on Slopes

*Bhargav Jyoti Borah and Sasanka Borah*

#### **Abstract**

Landslides in layered soil in slopes is a common phenomenon that leads to major damages to infrastructures. Piles are deep foundations which are useful structural elements to support heavy loads in difficult slopes, and to retain creeping or sliding slopes. Controlling the stability of structures made on slopes is still one of the major problems in engineering. Piles are being used successfully to support structures on slopes and have become an efficient solution for such situations which demands construction on sloping grounds, since piles can often be easily installed without disturbing the equilibrium of the layers of soil in the slope. In this paper, an optimized design approach have been adopted for design of single piles for supporting infrastructure on layered soil slopes. A standalone software has been developed as a result of this study which produces the optimum pile dimension, reinforcement details and an estimate of the design as software output. The code is written in MATLAB that incorporates Genetic Algorithm (GA) optimization technique to obtain optimum pile dimensions and the pile is designed satisfying all the structural requirements as recommended by IS:2911: Indian Standards Code for cast-in-situ Reinforced Concrete pile.

**Keywords:** piles on slopes, computer aided design of pile foundations, deep foundation, GA, MATLAB software, optimization algorithm

#### **1. Introduction**

Piles are deep foundations which are introduced into the soil by suitable means to support the load coming on it from the superstructure when a good bearing stratum is not available at shallow depths from the ground surface [1]. In such situations load have to be transmitted to a firm strata which is capable of supporting such tremendous loads at an appreciable deeper depth below the ground surface. Piled structures are often subjected to significant lateral loads in addition to the vertical load coming from the superstructure which are installed in soil layers present in sloping grounds. The varying properties of multi-layer soil in the slopes may further complicate the stability issue of the structure.

In this research, the model superstructure considered on multi-layer soil slope is as per Indian standards recommendation [2]. Indian standards guidelines have been followed for the design of pile foundation in the soil slope. Indian Standards [3–6]

classifies Concrete Piles as Driven Cast-In situ piles, Bored Cast In-situ piles, Driven Precast piles and Precast piles in Pre-bored Holes. The Load Bearing capacity of the pile is dependent on the material used and the dimensions of the pile foundation, the spacing of individual piles in the pile group, the load bearing capacity and nature of the supporting soil. The pile installation method, and the direction of incoming loading [7] also plays a significant role in determining the load bearing capacity of the pile foundation. The load bearing capacity of the soil slope is affected by many factors such as type of the soil in each layer of the slope and strength of soil in these soil layers, foundation dimensions, unit weight of the soil, surcharge acting on the soil, the type of loading etc. [1]. The variability of the properties of the supporting soil, soil profiles and properties of the soil in these multi-layered soil profiles, which generally exists in nature, affects the load bearing capacity of the pile drastically. The choice of a pile to be used is governed by the condition of the site, economy, time available etc. The problem of landslides is also a major consideration for stability of a structure on sloping grounds. Hence, in order to counter these problems which have a large number of variability associated with them, the aid of modern high-speed computers and software has become popular.

Conventional design method for pile foundations involves determining the load carrying capacity of pile foundations from Static method, Dynamic method, Pile load test method or Penetration Test method depending upon the situation giving due consideration to data obtained from soil exploration reports. The pile foundation or group of piles are designed in such a way so as to transfer the load to the supporting soil safely giving due consideration to the load coming from the superstructure. For the design consideration of geotechnical design of the pile, the Ultimate Load Carrying Capacity and the Allowable Load Carrying Capacity is taken into consideration whereas for the structural design of the foundation, the bearing capacity, driving and handling stresses of the pile foundation is taken into consideration. Design of Pile Foundation is generally a trial-and-error process. Here an initial trial design (foundation dimensions and reinforcements) is taken and is checked against the geotechnical and structural requirements. This is done iteratively to revise the trial designs till a practical design of pile foundation is obtained. The drawback of this general method is that with these methods, the design of the pile foundation often comes out to be a conservative one. Generally such a conservative design may lead to an uneconomical and an impractical design. In a general optimization approach for design of a pile foundation, an optimization algorithm is used. A well calibrated optimization method may confirm an economic and a practical pile foundation design. In optimized design method of design for single piles the Ultimate Bearing Capacity of the pile is evaluated within stipulated limits by changing the pile dimensions iteratively until an optimized pile dimension is found to support the total incoming load. The pile is then designed structurally for bearing the loads acting on the pile satisfying all the structural requirements, may it be for vertical loading or lateral loading.

Complex studies have been carried out based on experimental and numerical investigations of piles installed in sloping grounds. But the research done in foundation engineering for the application of optimization methods using Computer Aided Design methods for structures on sloping grounds is scarce. The available research that deals with pile design optimization, assume different methodology with a number of assumptions which cannot involve all the practical parameters and thus have limited practical applications for supporting structures on sloping grounds.

Borah and Borah [1] has adopted a GA based optimization algorithm for optimization and design of Vertical piles on horizontal grounds with the Ultimate bearing

*Genetic Algorithm Based Software for Optimization and Design of Piles on Slopes DOI: http://dx.doi.org/10.5772/intechopen.108615*

capacity being the objective function and the pile dimensions as the design variables. Hoback and Truman [8] have introduced a weightless optimality rule into the original optimality criteria approach to address the design variables- the spacing and battering of the piles that has no measurable effect on the objective function. Huang and Hinduja [9] adopted a quasi-Newton method for optimizing the shape of piles. Their algorithm includes a linear force–deflection relationship for the pile-soil system while optimizing the pile design. Nikolaou and Pitilakis [10] developed a stand-alone program based on MATLAB. The have included algorithm for the calculation of bearing capacity and settlements of shallow foundations using several well-known formulas. They have utilized various literature and codes of practice that are preferred in engineering.

In this paper, the research done in Borah and Borah [1] has been taken as the base for the development of a Genetic Algorithm based software code for optimization and design of concrete cast in situ pile in a multi-layer soil slope. The Genetic Algorithm (GA) is used to optimize the pile dimensions based on Indian Standards [4–7] methodology within recommended limits and user defined practical specifications. The piles designed by the software takes into account the lateral load and the lateral moment, along with the vertical loading and moment coming from the superstructure, which will be unique for piles on sloping grounds with multiple layers of soil. A software has been developed using MATLAB coding for automating the entire process of design and optimization. An estimate [11] for construction is also presented which may further ease the choice of selection of a particular design.

#### **2. Methodology**

#### **2.1 General**

In this study a genetic algorithm (GA) has been developed for optimization and design of a pile foundation on sloping ground. The Optimization method used for optimization in the said software is genetic algorithm (GA). Genetic algorithm (GA) is a search-based optimization technique that works on the principles of Genetics and Natural Selection which is used to find optimal solutions to difficult problems which otherwise would take a considerable amount of time to solve [1]. It is generally used for solving optimization problems in research. The Pile Design, both geotechnical and structural, is as per Indian standards code [4–7, 12, 13] based design method. A preliminary estimate of conceptual design is also presented which is as per Assam Public Works Department Schedule of Rates [11]. MATLAB which is a programming software has been used to develop the interface of the stand-alone software and codes in this study.

#### **2.2 Genetic algorithm**

The Genetic Algorithm (GA) is a stochastic search-based optimization algorithm which mimics the process of evolution that includes natural selection. It is commonly used for nonlinear optimization problems. On the other hand, gradient based methods can be applied to linear optimizations due to certain defects in accuracy of the objective function [1]. The Genetic algorithm involves processes of selection, crossover, and mutation of the trial solutions so as to improve the set of solutions and converge towards an optimal solution. Each solution in every generation consists of a set of

parameters, and by manipulating these parameters the Genetic Algorithm (GA) converges towards the best possible solution. In Genetic Algorithm (GA), a population of solutions, also referred to as individuals are generated. For a given optimization problem, a set of potential solutions are initially generated, also known as the first generation. Within each population, every solution or individual is assigned a corresponding fitness value which is obtained by calculating a fitness function. The fitness value ranks the fitness or appropriateness of a solution. Within each population, some solutions based on their fitness are chosen to carry forward to the next generation of solutions as clones of the original solution. These solutions thus guarantee with the best fitness value of each generation will be maintained or manipulated to be improved upon from one generation to another, while providing better quality of said parameters within a population.

Crossover represents a reproduction function. It means that for each generation, a smaller group of solutions is selected to combine the parameters associated with them and create new solutions, and these new solutions are the new generation of solutions. The group that crosses over are called parents, and the "genes" of two (or more) such parent solutions are combined to generate one (or more) new solution, as their child. It is these children that make up a new generation of solutions. These new generation of solution is expected to be better in terms of fitness as they are made from parents having best fitness in their respective generation [1].

Each solution are then subjected to mutations. It is a random change in one or more parameters (genes) of a solution based on some probability. This probability is known as the mutation rate [1]. Mutations is a must to maintain and introduce diversity into a population of solution. With more diversity in the population, the genetic algorithm lowers the risk of ending up with a solution in the sub-optimal local minima. A flowchart of the simple version of Genetic Algorithm (GA), explaining its process is shown in **Figure 1**.

#### **2.3 Design considerations**

IS 2911(Part 1/Sec 2): 2010 [5] recommends that the design of a pile foundation should be in such a way that the incoming load from the super structure can be transmitted to the sub-surface safely with appropriate factor of safety. This Factor of safety is against shear failure of sub-surface and without causing excessive settlement (differential or total). The shaft of the pile should have adequate structural capacity to withstand all types of incoming loads, may it be vertical, axial or otherwise and incoming moments which is to be transmitted to the supporting subsoil. The Pile should be structurally complying with the design recommendations given in IS 456: 2000 [13].

#### **2.4 Pile capacity**

In this study the ultimate load capacity of a pile foundation is obtained by using a static analysis methodology recommended by IS 2911(Part 1/Sec 2): 2010 [5]. The Load Carrying Capacity of the pile depends on the of the supporting soil properties of various layers of the supporting soil in which the pile is installed. The minimum factor of safety on static formula is taken to be 2.5 as is recommended by the code, IS 2911 (Part 1/Sec 2): 2010 [5].

*Genetic Algorithm Based Software for Optimization and Design of Piles on Slopes DOI: http://dx.doi.org/10.5772/intechopen.108615*

**Figure 1.** *Basic genetic algorithm (GA) flow chart.*

#### **2.5 Analysis of laterally loaded piles**

A pile foundation which is subjected to lateral forces from a number of sources, such as, wind, earthquake, water current, earth pressure, effect of moving vehicles or ships, plant and equipment, etc. [5]. The lateral load carrying capacity of a pile foundation depends on the horizontal sub-grade modulus of the supporting soil and on the structural capacity of the shaft of the pile against bending. While considering lateral load on pile foundations, the effect of other loads acting on the pile foundation, like the axial load coming from the super structure, is to be taken into consideration.

The IS code [5] suggests that a group of three or more pile connected by a rigid pile cap is considered to have a fixed head condition. In all other cases the pile foundation is taken to be free headed.

#### **2.6 Structural capacity**

The IS code [5] suggests that the pile foundation should have necessary structural strength to transmit the imposed loads safely and ultimately to the supporting soil.

#### **2.7 Estimation**

An estimate of the pile foundation designed by the algorithm is generated by the software. In this study the software generates the design of the pile foundation as the output which includes the structural details of the pile foundation and various specifications that will be required in its construction. The structural design generated by the software serves as the drawing for estimation, the other concrete and steel parameters along with workmanship requirements specified by Indian Standards [5] are considered and the Schedule of Rates [11] provided by the Assam Public Works Department, provides the Rate of Items associated with the design.

#### **3. Results and discussions**

#### **3.1 Software**

The software which is developed in this research has been designed to optimize pile dimensions based on the incoming load (axial, eccentric and lateral) from the super structure and the supporting soil parameters of sloping ground. The software developed allows user to use optimized results generated by the software to design the pile foundation structurally or alternatively allow user input data for design of the pile foundation. All pile foundation design procedure has been adopted as per guidelines provided by Indian Standards code [4–7, 12, 13] on sloping grounds. A detailed step by step procedure for using this software for optimization and design of pile foundation based on the parameters of soil layers on sloping ground is given below:


*Genetic Algorithm Based Software for Optimization and Design of Piles on Slopes DOI: http://dx.doi.org/10.5772/intechopen.108615*



A stand-alone software is coded with the use of MATLAB Compiler to allow users without a MATLAB license to install and run the software on different operating platforms.

The entire process is summarized in the flow diagram (see **Figure 3**).

#### **3.2 Calibration of software**

A few numerical examples [14–16] were taken from for calibrating the software developed in this study. Methodology discussed in the previous section has been adopted for optimization and design of the pile foundation on sloping ground. The results obtained by the software, by using the parameters given in these numerical problems, are then verified manually.

The basic form of the Genetic Algorithm (GA) used while developing the software is as follows:

**Optimize:**

$$\mathbf{Q}\_{\mathbf{u}} = \mathbf{A}\_{\mathbf{p}} \left( \frac{\mathbf{1}}{2} \mathbf{D}\_{\mathbf{y}} \mathbf{N}\_{\mathbf{y}} + \mathbf{P}\_{\mathbf{D}} \mathbf{N}\_{\mathbf{q}} \right) + \sum\_{i=0}^{n} \mathbf{K}\_{i} \mathbf{P}\_{\mathbf{Di}} \tan \delta\_{i} \mathbf{A}\_{\mathbf{si}} \tag{1}$$

$$\text{for } \mathbf{Q}\_{\text{u}} = \mathbf{A}\_{\text{p}} \mathbf{N}\_{\text{C}} \mathbf{c}\_{\text{p}} + \sum\_{\text{i=0}}^{\text{n}} \alpha\_{\text{i}} \mathbf{c}\_{\text{i}} \mathbf{A}\_{\text{si}} \tag{2}$$

Based on the nature of the soil layers on sloping ground. The first term is the endbearing resistance (Qp) and the second term is the skin friction resistance (Qs).

#### **Subject to the constraints:**


*Genetic Algorithm Based Software for Optimization and Design of Piles on Slopes DOI: http://dx.doi.org/10.5772/intechopen.108615*

#### **Figure 3.**

*Flow diagram depicting pile foundation optimization and design process.*


The results of software developed in this study are verified manually and the results tallies. Hence the working process of the software has been calibrated with those of a standard.

#### **3.3 Findings and interpretations**

Based on the Results obtained from the software in this study and subsequent manual verification of the results, it can be inferred that (**Figures 4**–**6**):

#### **Figure 4.**

*Optimization and design results for three layers of soil.*

#### **Figure 5.**

*Optimization and design results for two layers of soil.*


*Genetic Algorithm Based Software for Optimization and Design of Piles on Slopes DOI: http://dx.doi.org/10.5772/intechopen.108615*

**Figure 6.** *Optimization and design results for single layer of soil.*


#### **4. Conclusion**

The software developed through this study has adopted the design methodology and constraints considering all the variables involved in pile foundation design and construction recommended by Indian Standards code [4, 5, 6, 7, 12], which is widely used in practice in India. The software developed in this study lets user generate designs for an optimized pile foundation on sloping grounds. The software can also be used for generating structural design of a pile foundation subjected to lateral loading in addition to vertical loading. The said software takes into account the variability of soil layer parameters on sloping grounds thus enhancing the scope for further research into the code for countering complex practical problems. Due to its simplicity and user friendly interface it may cater to shorten the gap between conceptual optimization and design outputs and the actual field applicability of piles on slopes. The software also provides an estimate of the design produced so as to allow the designer to better judge the practical applicability of a particular design. Through this study an attempt has been made to bridge the gap between research and field application of an optimized foundation design on sloping grounds and incorporating computer coding to develop an application to simplify user experience.

*Current Perspectives on Applied Geomorphology*

### **Author details**

Bhargav Jyoti Borah<sup>1</sup> \* and Sasanka Borah<sup>2</sup>

1 Chirang Polytechnic, Chirang, India

2 Assam Engineering College, Guwahati, India

\*Address all correspondence to: borahjbhargav@gmail.com

© 2023 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

*Genetic Algorithm Based Software for Optimization and Design of Piles on Slopes DOI: http://dx.doi.org/10.5772/intechopen.108615*

#### **References**

[1] Borah BJ, Borah S. Genetic algorithm based optimization and Design of Pile Foundation. In: Satyanarayana Reddy CNV, Muthukkumaran K, Satyam N, Vaidya R, editors. Ground Characterization and Foundations. Lecture Notes in Civil Engineering. Vol. 167. Singapore: Springer; 2022

[2] IS 14343 (Part 2): Seclection and Development of Site for Building in Hill Areas — Guidelines —Selection and Development. New Delhi: Bureau of Indian Standards; 1995

[3] IS 2911 (Part 1/Sec 1): Design and Construction of Pile Foundations — Code of Practice Part 1 Concrete Piles Section 1 Driven Cast In-Situ Concrete Piles. New Delhi: Bureau of Indian Standards; 2010

[4] IS 2911 (Part 1/Sec 2): Design and Construction of Pile Foundations — Code of Practice Part 1 Concrete Piles Section 2 Bored Cast In-Situ Concrete Piles. New Delhi: Bureau of Indian Standards; 2010

[5] IS 2911 (Part 1/Sec 3): Design and Construction of Pile Foundations — Code of Practice Part 1 Concrete Piles Section 3 Driven Precast Concrete Piles. New Delhi: Bureau of Indian Standards; 2010

[6] IS 2911 (Part 1/Sec 4): Design and Construction of Pile Foundations — Code of Practice Part 1 Concrete Piles Section 4 Precast Concrete Piles in Pre-Bored Holes. New Delhi: Bureau of Indian Standards; 2010

[7] Abdurrahman S. Mathematical models and solution algorithms for computational design of RC piles under structural effects. Applied Mathematical Modelling. 2011;**35**(1):3611-3638

[8] Hoback AS, Truman KZ. Least weight design of steel pile foundations. Engineering Structures. 1993;**15**(5): 379-385

[9] Huang Z, Hinduja S. Shape optimization of a foundation for a large machine tool. International Journal of Machine Tool Design and Research. 1986;**26**(2):85-97

[10] Nikolaou K, Pitilakis D. SoFA: A matlab-based educational software for the shallow foundation analysis and design. Computer Applications in Engineering Education. 2017;**25**(2): 214-221

[11] Assam Public Works Department Schedule of Rates 2013–14 (Building). Assam: Commissioner & Special Secretary, Public Works Department; 2013

[12] IS 6403: Indian Standard Code of Practice for Determination of Breaking Capacity of Shallow Foundations. New Delhi: Bureau of Indian Standards; 1981

[13] IS 456: Plain and Reinforced Concrete Code of Practice. New Delhi: Bureau of Indian Standards; 2000

[14] Murthy VNS. Textbook of Soil Mechanics and Foundation Engineering. 2nd ed. New Delhi: CBS Publishers & Distributors; 2009

[15] Punmia BC, Jain AK, Jain AK. Soil Mechanics and Foundations. 16th ed. Guwahati: Laxmi Publications (P) Ltd.; 2014

[16] Ranjan G, Rao ASR. Basic and Applied Soil Mechanics. 2nd ed. Guwahati: New Age International Publishers (P) Ltd.; 2014

[17] SP:16: Design Aids for Reinforced Concrete to IS: 456-l978. New Delhi: Bureau of Indian Standards; 1987

[18] SP:34 (S&T): Handbook on Concrete Reinforcement and Detailing. New Delhi: Bureau of Indian Standards; 1987

#### **Chapter 3**

### Slope Stability Analyses Subjected to Slide Head Toppling Failure Mechanisms

*Victor Mwango Bowa and Abdul Samson*

#### **Abstract**

It is faulty to analyze the jointed rock slopes'stability susceptible to a combination of modes of failure composed of sliding around the toe region and toppling of the rock blocks on the upper part of the slope based on the current analytical methods, which are based on assumption that the distribution of the potential failure surface bounding the potential mixed failure runs predictably from crest to toe of the slope. An Analytical model that takes into consideration the kinematic mechanism of the discontinuous rock slope with counter-tilted weak plane subjected to a combination of failure mechanisms involving sliding and toppling has hence been presented in this chapter. This involves an iterative process which involves the calculations of the dimensions of all the individual blocks as well as the forces acting on them, and then stability of every block is examined, starting at the uppermost block. The stability analysis of each block is determined. The blocks may either be stable, topple or slide. The proposed analytical methods could curtail errors incurred due to the acceptance of the single weak plane for quantifying the failure mechanisms composed of slide head toppling rock slopes in physical situations with two planar weak planes.

**Keywords:** failure mechanisms, jointed rock slopes, stabilizing techniques, analytical solution, slide head toppling

#### **1. Introduction**

A lot of research has been done in the current analytical models for predicting rock slope stability subjected to failure by block toppling. However, their contributions have focused always on the predictive and idealized geometry comprising blocks with joints dipping into the slope face. In addition, the jointing bounding the latent toppling blocks is predicted to be systematic running from the uppermost surface of the slope, and day-lighting near the toe. It is however possible that in some notable situations, the jointing (discontinuity) bounding the potential toppling rock blocks may not be predictive due to geological variations. Hence, the joints that bound the potential sliding and toppling blocks may be counter-tilted in the rock mass due to geological variations and may daylight on an unpredicted point on the

slope face. This would lead to a combination of failure mechanisms composed of toppling and other secondary failure mechanisms such as toppling-circular, slidehead-toppling and block-flexure toppling. Not much research has been done on a combination of failure modes which involve sliding and block toppling on the discontinuous rock slope.

Regardless of this, a combination of failure mechanisms involving sliding and block toppling of the upper section of the slope respectively, continue to be a common phenomenon in sedimentary rock formation in numerous constructed and natural slopes all over the world as illustrated in **Figure 1a,** and **b**.

Furthermore, another typical example of slide-head-toppling failure is an open pit mine slope in Valencia in Spain shown in **Figure 2**. **Figure 2** shows a series of stable rock blocks in the soil mass on the upper surface of the slope, another set of toppling columnar blocks in the limestone formation midway section and finally a series of sliding blocks in the variation of the weak rock formation on the slope's lower part.

It is not entirely correct to analyze the stability of discontinuous rock slopes with a high potential to mixed toppling failure modes of toppling on the upper section of the slope and sliding on the slope's lower section using the existing limit equilibrium method, which is based on the initial assumptions of distribution of the potential failure surface bounding the blocks susceptible mixed failure mode. The theoretical/analytical model is therefore necessary that it takes into consideration the kinematic mechanisms of the discontinuous rock slopes composed of the counter-tilted jointing susceptible/subjected to a combination of failure mechanisms (sliding and toppling).

This chapter presents a two-dimensional (2D) model which incorporates the counter-tilted oriented jointing bounding the potential for both sliding and toppling. From the 2D model, a theoretical (analytical) model that is based on the trigonometry basic principles to determine the geometry of the jointed rock slope composed of the

#### **Figure 1.**

*Mixed toppling failure modes involving toppling and sliding failures observed in practice (a) post-failure photograph of sliding and toppling of north wall slope of Nchanga open pits, Zambia [1] (b) post-failure photograph of sliding and toppling failure in Zhongliang reservoir bank, China [2].*

*Slope Stability Analyses Subjected to Slide Head Toppling Failure Mechanisms DOI: http://dx.doi.org/10.5772/intechopen.108762*

**Figure 2.**

*Valencia open pit mine with the potential of mixed failure modes composed of toppling and sliding failures in Spain [3].*

counter-tilted jointing subjected to slide-head toppling failure mechanisms have been put forward. Next, a process of iteration is followed with forces acting on each block determined. A comparative analysis between the stability against sliding and toppling of individual blocks is determined, with the highest forces between these two deciding the failure mechanism of the succeeding block. In summary, the overall slope stability is considered based on whether the lowermost block either slides or topples.

#### **2. Rock slope geometry**

For stability analyses of the slope to be achieved, detailed slope geometries must be determined at the initial stage. Let us examine a rock slope in **Figure 3**, which is composed of rectangular rock columns with width, Δ*x* and height, *yn* in an orderly manner. The original failure plane dipping at *ψp* is counter-tilted from the preliminary weak plane approximately on the mid-way part of the slope at angle, *θr*. The counter-tilted failure plane then dips at *ψc*. Other slope parameters include the slope height indicated as H, the face angle denoted as *ψf* and the dip/angle of the uppermost slope face referred to as *ψs*. The rock block column count started from the lower most of the slope (toe) increasing cumulatively upwards. It has been observed and noted from centrifugal and numerical test models that the base plane tends to be stepped during toppling failure mechanisms [4–8]. This greatly influences the determination of the overall angle of the plane base denoted as *ψb*.

In rock slopes where the failure plane is counter-tilted around the midway section of the rock slope, it is important to approximate the base plane angle (*ψb*) considering

**Figure 3.** *Rock slope composed of counter-tilted weak surface subjected to toppling failure mechanism [1, 4].*

the dipping of the failure plane bases for the rock block columns at *ψp* and *ψc* using Eq. (1)

$$
\mu \nu b \approx \left(\frac{\nu c + \nu p}{2} + 10^{\circ}\right) \text{ Or} \nu b \approx \left(\frac{\nu c + \nu p}{2} + 30^{\circ}\right) \tag{1}
$$

Based on the slope geometry in **Figure 3**. Eq. (2) is formulated to determine the number of rock columns, n, forming the regular system of the slope.

$$n = \frac{H}{\Delta x} \left[ \cos \alpha c(yb) + \left( \frac{\cot \left( yb \right) - \cot \left( yf \right)}{\sin \left( yb - yf \right)} \right) \sin \psi s \right] \tag{2}$$

Considering the slope model presented in **Figure 3**, the height, *yn* for the nth rock column below the slope's crest dipping/orientated into the slope face is determined by Eq. (3).

$$yn = n(\mathbf{a1} - b) \tag{3}$$

The height, *yn* for the rock column denoted as nth above the slope's crest dipping into the slope face is resolved by Eq. (4).

$$
\gamma m = \gamma n - \mathbf{1} - \mathbf{a}\mathbf{2} - b \tag{4}
$$

The constants a1, a2 and *b* shown in **Figure 3** are found based on the rock columns and the related slope geometry. These constants are calculated by Eqs. (5)–(7) respectively;

*Slope Stability Analyses Subjected to Slide Head Toppling Failure Mechanisms DOI: http://dx.doi.org/10.5772/intechopen.108762*

$$\mathbf{a1} = \Delta \mathbf{x} \tan \left( yf - \psi c \right) \tag{5}$$

$$\mathbf{a}\mathbf{2} = \Delta\mathbf{x}\tan\left(\wp p - \wp s\right) \tag{6}$$

$$b = \Delta \mathfrak{x} \tan\left(\psi b - \psi p\right) \tag{7}$$

Other distinctive dimensions of the blocks required to be determined includes; the points of application denoted *Mn* and *Ln* for the shear and normal forces (Rn,Sn) that acts on the bases of weak planes and (Pn,Qn,Pn � 1,Qn � 1) forces applied on interfaces next to the rock columns during toppling and sliding failure mechanisms as illustrated in **Figure 3**. The application points of all forces needs to be determined every time toppling failure mechanisms takes place. Considering the case where the nth block is below the slope's crest, the application points *Mn* and *Ln* are determined by Eqs. (8) and (9) respectively, at the crest block, then the application points *Mn* and *Ln* are determined by Eqs. (10) and (11) and finally when the nth block is over the slope's crest, then the application points *Mn* and *Ln* are calculated by Eqs. (12) and (13).

$$Mn = \text{yn} \tag{8}$$

$$Ln = \mathcal{ym} - \mathbf{a1} \tag{9}$$

#### **Figure 4.**

*Limit equilibrium conditions required for modes of failure of the nth rock columns: (a) forces applied on the nth rock columns; (b) toppling mode of failure of the nth rock columns; (c) sliding modes of failure of the nth rock columns considering the counter-tilted failure plane.*

$$Mn = \text{yn} - \text{a}\mathbf{2} \tag{10}$$

$$Ln = \mathfrak{yn} - \mathfrak{a}\mathbf{1} \tag{11}$$

*Mn* ¼ *yn* � a2 (12)

$$Ln = \mathfrak{yn} \tag{13}$$

#### **2.1 Block stability**

**Figure 4** illustrates a rock slope composed of the counter-tilted weak (failure) surface composed of three sets of rock blocks that have been categorized based on their failure behavior. The rock columns for the upper most section of the slope are regarded as stable based on the fact that the rock columns' base friction angle is more than the original failure plane angle. For the rock block columns in the midway section of the slope, toppling failure is observed based on the consideration that the base plane lies outside the Centre of gravity for the rock block columns. For the blocks located around the lowermost section of the slope, there is a high possibility of rock blocks undergoing failure as the initial failure plane is counter-tilted. The modes of failure for the block columns are influenced by the geometry of the slope, the slope angle and the failure/weak plane angle. The friction angles of these failure surfaces differ with respect to characteristics of the lithology.

Generally, for such conditions, the friction angle of the block interfaces' of the rock block columns denoted (*ϕd*) are presumed to be equal to the angle of friction on the base plane (*ϕ*p,*ϕc*). The two sides'shearing forces of the rock column are then resolved using Eqs. (14) and (15), considering the limiting equilibrium analysis. On the other hand, the forces applied on the rock column considering the orthogonal and sides' friction angles to the two base planes (*ψp*, *ψc*) of the rock column having weight, *Wn*, are resolved using Eqs. (16) and (17).

$$Qn = Pn \tan \phi d \tag{14}$$

$$Qn - 1 = Pn - 1 \tan \phi d \tag{15}$$

$$\begin{aligned} Rn &= Wn\cos\psi p + (Pn - Pn - 1)\tan\phi d \\ Rn &= Wn\cos\psi c + (Pn - Pn - 1)\tan\phi d \end{aligned} \tag{16}$$

$$\begin{aligned} \text{Sn} &= \text{Wn}\sin\psi p + (Pn - Pn - \mathbf{1})\\ \text{Sn} &= \text{Wn}\sin\psi c + (Pn - Pn - \mathbf{1}) \end{aligned} \tag{17}$$

To define the magnitude of *Pn*, the moment force about the pivot point O is set to zero (refer **Figure 4**). We then, as discussed above, assume that the angle of friction of the interfaces of the rock block columns denoted (*ϕd*) is equal in magnitude to the angle of friction on the base planes (*ϕ*p, *ϕc*), thus *ϕd* ¼ *ϕp* ¼ *ϕc* ¼ *ϕ*.

Considering the possibilities of rotational equilibrium, the forces Pn,t necessary to avoid toppling of the *nth* block considering the presumed angle of the failure plane is resolved based on Eq. (18) considering the derivations as follows:

The moment about the pivot point 0 is set to zero. Then we have:

$$
\sum \mathbf{M} \mathbf{0} \boldsymbol{\oplus} = \mathbf{0}
$$

*Slope Stability Analyses Subjected to Slide Head Toppling Failure Mechanisms DOI: http://dx.doi.org/10.5772/intechopen.108762*

$$\rho = \frac{\eta n}{2} W n \sin \psi p - \frac{\Delta \chi}{2} W n \cos \psi p + M n P n + 1 - \Delta \chi (P n + \mathbf{1} \tan \phi) - L n P n$$

$$= \frac{W n}{2} (\eta n \sin \psi p - \Delta \chi \cos \psi p) + P n + \mathbf{1} (M n - \Delta \chi \tan \phi) - L n P n$$

We can rewrite this as follows:

$$LnPn = \frac{Wn}{2} \left( yn \sin \psi p - \Delta x \cos \psi p \right) + Pn + \mathbf{1} (Mn - \Delta x \tan \phi)$$

Hence this gives the following equation for *Pn*

$$Pn, t = \frac{\frac{\eta\_n}{2} (yn \sin \psi p - \Delta \varkappa \cos \psi p) + Pn + \mathbf{1} (Mn - \Delta \varkappa \tan \phi)}{Ln} \tag{18}$$

Considering the possibilities of rotational equilibrium, the forces Pn,t necessary to avoid toppling for the *nth* block taking into consideration the counter-tilted angle of the failure plane is determined using Eq. (19) considering the derivations as follows;

Again, the moment about the pivot point O is set to zero. Then we have:

$$\sum M0\Theta = 0$$

$$\eta = \frac{\eta n}{2} Wn \sin \psi c - \frac{\Delta \chi}{2} Wn \cos \psi c + MnPn + \mathbf{1} - \Delta \chi (Pn + \mathbf{1} \tan \phi) - Ln Pn$$

$$= \frac{Wn}{2} \left( \eta n \sin \psi c - \Delta \chi \cos \psi c \right) + Pn + \mathbf{1} (Mn - \Delta \chi \tan \phi) - Ln Pn$$

We can rewritten this as follows;

$$LnPn = \frac{Wn}{2} \left( \mathcal{ym} \sin \psi c - \Delta \mathfrak{x} \cos \psi c \right) + Pn + \mathbf{1} (Mn - \Delta \mathfrak{x} \tan \phi)$$

Therefore, we obtain the following equation for *Pn*,*t* with respect to counter-tilting of the failure plane:

$$Pn, t = \frac{\frac{\mathcal{W}n}{2}(yn\sin\psi c - \Delta x \cos\psi c) + Pn + \mathbf{1}(Mn - \Delta x \tan\phi)}{Ln} \tag{19}$$

If the rock columns' failure mode on the slope's lower section shown in **Figure 4c** is noted to be sliding, then to calculate the magnitude of Pn,s, the total (sum) of the forces in the horizontal directions and vertical directions is set to zero (assuming the horizontal axis to lie along the base plane, at angles *ψp* or *ψc* considering the horizontal as the datum refer **Figure 4**. In general, it is presumed that the interfaces' angle of friction for the rock columns denoted (*ϕd*) is equal to the friction angle on the base plane (*ϕ*p, *ϕc*), thus *ϕd* ¼ *ϕp* ¼ *ϕc* ¼ *ϕ*.

Considering the possibilities of rotational equilibrium, the forces Pn,*s* required to prevent sliding of block n considering the initial angle of the weak plane is resolved using Eq. (20) considering the derivations as follows;

Before counter-tilting of the failure plane along base plane angle *ψp* taking the horizontal as the datum, we have

$$
\sum F\mathbf{x} = \mathbf{0} = Pn - Pn + \mathbf{1} - Wn\sin\psi p + Sn
$$

$$
\sum F\mathbf{y} = \mathbf{0} = Rn + Pn\tan\phi - Wn\cos\psi p - Pn + \mathbf{1}\tan\phi
$$

Where *Sn* is shear force acting along the base of the column contacts, from the Mohr-Coulomb criterion, we have:

$$\pi a = c + \sigma a \tan \phi$$

Where *τn* shear stress acting along the base of the column contact and *σn* Normal stress at the base column contact.

We further assumed the cohesion (c) being negligible and divide both sides of the equation by the column-base plane contact area denoted by A*.* Thus, we get

$$\mathfrak{Sn} = Rn \tan \phi$$

Where; *Rn* is normal force acting across the base plane of the column contacts (*Wn* cos *ψp*).

By resolving the *Fy* equation for *Rn* making substitutions into the *Fx* equation prior to counter-tilting of the failure plane along base plane angle *ψp* with regards to the horizontal, we have,

$$(Pn - Pn + \mathbf{1} - Wn\sin\psi p + Wn\cos\psi p\tan\phi - (Pn - Pn + \mathbf{1})\tan^2\phi = \mathbf{0})$$

$$(Pn - Pn + \mathbf{1})(\mathbf{1} - \tan^2\phi) - Wn\sin\psi p + Wn\cos\psi p\tan\phi = \mathbf{0}$$

This leads to:

$$(Pn - Pn + 1)(1 - \tan^2 \phi) = Wu \cos \psi p (\tan \psi p - \tan \phi)$$

$$(Pn - Pn + 1) = \frac{Wn \cos \psi p (\tan \psi p - \tan \phi)}{1 - \tan^2 \phi}$$

Hence, we obtain the following equation for *Pn*,*s*

$$Pn\,\varphi = Pn + \mathbf{1} + \frac{\mathcal{W}n\,\cos\varphi p(\tan\varphi p - \tan\phi)}{\mathbf{1} - \tan^2\phi} \tag{20}$$

Considering the possibilities of rotational equilibrium, the forces Pn,*s* necessary to avoid sliding of *nth* block taking into consideration the counter-tilted angle of the weak plane is resolved using Eq. (21) considering the derivations as follows;

Once counter-tilting of the failure plane takes place along base plane *ψc* taking into consideration the horizontal line, we have

$$
\sum F\mathbf{x} = \mathbf{0} = Pn - Pn + \mathbf{1} - Wu\sin\varphi c + Sn
$$

$$
\sum F\mathbf{y} = \mathbf{0} = Rn + Pn\tan\phi - Wu\cos\varphi c - Pn + \mathbf{1}\tan\phi
$$

By solving the *Fy* equation for *Rn* and making substitutions into the *Fx* equation once counter-tilting of the failure plane takes place along base plane *ψc* taking the horizontal as the datum, we have,

*Slope Stability Analyses Subjected to Slide Head Toppling Failure Mechanisms DOI: http://dx.doi.org/10.5772/intechopen.108762*

$$(Pn - Pn + \mathbf{1} - Wn\sin\psi c + Wn\cos\psi c\tan\phi - (Pn - Pn + \mathbf{1})\tan^2\phi = \mathbf{0})$$

$$(Pn - Pn + \mathbf{1})\left(\mathbf{1} - \tan^2\phi\right) - Wn\sin\psi c + Wn\cos\psi c\tan\phi = \mathbf{0}$$

This leads:

$$(Pn - Pn + \mathbf{1}) \left(\mathbf{1} - \tan^2 \phi\right) = Wn \cos \psi c (\tan \psi c - \tan \phi)$$

$$(Pn - Pn + \mathbf{1}) = \frac{Wn \cos \psi c (\tan \psi c - \tan \phi)}{\mathbf{1} - \tan^2 \phi}$$

Hence, the following equation for *Pns* is obtained;

$$Pn\varphi = Pn + \mathbf{1} + \frac{Wn\cos\varphi c(\tan\varphi c - \tan\phi)}{\mathbf{1} - \tan^2\phi} \tag{21}$$

#### **2.2 Determination of anchor tension necessary to prevent toppling of block 1**

Once a determination is made that block 1 in **Figure 4c** will topple then cables under tension can be installed through block 1 to be anchored in stable rockmass underneath the toppling zone to prevent toppling of the same. Further assumptions was made that the anchor is installed at an angle plunge of *ψT* through block 1 located at a distance *L*1 above the base. After the application of force *T* to block 1, the normal force ð Þ *R*1 as well as shear force ð Þ *S*1 on the base plane of the block 1 is determined using Eqs. (22) and (23) respectively.

$$R1 = P1\tan\phi + T\sin\left(\psi p + \psi T\right) + W1\cos\psi p \tag{22}$$

$$\mathbf{S1} = P\mathbf{1} - T\cos\left(\varphi p + \varphi T\right) + W\mathbf{1}\sin\varphi p \tag{23}$$

To solve for the magnitude of *Tt*, again the moment force about the pivot point 0 is set to zero as follows:

$$\sum \mathbf{M} \mathbf{0} \oplus \mathbf{0}$$

$$\mathbf{0} = \frac{W \mathbf{1}}{2} (y \mathbf{1} \sin \psi p - \Delta x \cos \psi p) + P \mathbf{1} (y \mathbf{1} - \Delta x \tan \phi) - T t L \mathbf{1} \cos \left(\psi p + \psi T\right))$$

This is resolved as follows:

$$\frac{TtL1\cos\left(\wp p + \wp T\right)}{L1\cos\left(\wp p + \wp T\right)} = \frac{\frac{W1}{2}\left(\mathfrak{y}\mathbf{1}\sin\wp p - \Delta\mathfrak{x}\cos\wp p\right) + P\mathbf{1}\left(\mathfrak{y}\mathbf{1} - \Delta\mathfrak{x}\tan\phi\right)}{L1\cos\left(\wp p + \wp T\right)}$$

The magnitude of *Tt*, is deduced as shown in Eq. (24)

$$Tt = \frac{\frac{W1}{2}(y\mathbf{1}\sin\psi p - \Delta x\cos\psi p) + P\mathbf{1}(y\mathbf{1} - \Delta x\tan\phi)}{L\mathbf{1}\cos\left(\psi p + \psi T\right)}\tag{24}$$

However, if there is existence of the counter-tilted failure plane at angle *ψc* to the horizontal plane, the normal force ð Þ *R*1 and shear force ð Þ *S*1 on the base plane of the block 1 are resolved using Eqs. (25) and (26) respectively.

$$R1 = P1\tan\phi + T\sin\left(\psi c + \psi T\right) + W1\cos\psi p \tag{25}$$

$$\mathbf{S1} = P\mathbf{1} - T\cos\left(\psi c + \psi T\right) + W\mathbf{1}\sin\psi c \tag{26}$$

To solve for the magnitude of *Tt*, again the moment about the pivot point 0 is set to zero;

$$\mathbf{0} = \frac{W\mathbf{1}}{2}(\mathbf{y}\mathbf{1}\sin\psi\mathbf{c} - \Delta\mathbf{x}\cos\psi\mathbf{c}) + P\mathbf{1}(\mathbf{y}\mathbf{1} - \Delta\mathbf{x}\tan\phi) - T\mathbf{t}L\mathbf{1}\cos\left(\psi\mathbf{c} + \psi T\right)\mathbf{1}$$

This is resolved this as follows:

$$\frac{TtL1\cos\left(\wp c + \wp T\right)}{L1\cos\left(\wp c + \wp T\right)} = \frac{\frac{W1}{2}\left(\wp \mathbf{1}\sin\wp c - \Delta\mathbf{x}\cos\wp c\right) + P\mathbf{1}\left(\wp \mathbf{1} - \Delta\mathbf{x}\tan\phi\right)}{L1\cos\left(\wp c + \wp T\right)}$$

The magnitude of *Tt*, is reduced to Eq. (27).

$$Tt = \frac{\frac{W1}{2}(y\mathbf{1}\sin\psi c - \Delta x\cos\psi c) + P\mathbf{1}(y\mathbf{1} - \Delta x\tan\phi)}{L\mathbf{1}\cos\left(\psi c + \psi T\right)}\tag{27}$$

#### **2.3 Anchor tension necessary to prevent sliding of block 1**

After a determination is made of the sliding block 1in **Figure 4c** the cables under tension can therefore be fixed through block 1 and anchored in the stable rock mass underneath the zone of sliding to avoid movements of block 1. It is also assumed that the installed anchor at a plunge angle of *ψT* through block 1 at a distance *L*1 above the base plane. Therefore, with the failure plane along base angle *ψp* taking the horizontal as the datum, P*Fx* and P*Fy* are resolved using Eqs. (28) and (29) respectively.

$$
\sum F \mathbf{x} = \mathbf{0} = Pn - Pn + \mathbf{1} - Wu \sin \varphi p + \mathbf{S}n\tag{28}
$$

$$\sum Fy = 0 = Rn + Pn\tan\phi - Wn\cos\psi p - Pn + \mathbf{1}\tan\phi\tag{29}$$

Once the force *T* is applied to block 1, the normal force ð Þ *R*1 and shear force ð Þ *S*1 on the base planer of the block 1 are determined using Eqs. (30) and (31) respectively.

$$R1 = P1 \tan \phi + T \sin \left(\varphi p + \varphi T\right) + W1 \cos \varphi p \tag{30}$$

$$\mathbf{S1} = P\mathbf{1} - T\cos\left(\wp p + \wp T\right) + W\mathbf{1}\sin\wp p \tag{31}$$

By solving for the *Fy* equation for *Rn* and making substitution into the *Fx* equation prior to counter-tilting of the failure plane along the base plane *ψp* taking the horizontal line as the reference point, hence, the moments about the pivot point 0 was set to zero as given below;

$$\begin{aligned} 0 &= P\mathbf{1}(\mathbf{1} - \tan\phi p \tan\phi - W\mathbf{1}(\tan\phi p \cos\psi p - \sin\psi p) - Ts(\tan\phi \sin\psi p + \psi T) \\ &- \cos\left(\psi p + \psi T\right) \end{aligned}$$

This was resolved as follows:

$$\operatorname{Ts}(\tan\phi\sin\left(\psi p + \psi T\right) - \cos\left(\psi p + \psi T\right) = P\mathbf{1}(1 - \tan\phi p \tan\phi - W\mathbf{1}(\tan\phi p \cos\psi p - \sin\psi p))$$

*Slope Stability Analyses Subjected to Slide Head Toppling Failure Mechanisms DOI: http://dx.doi.org/10.5772/intechopen.108762*

$$\frac{\operatorname{Ts}(\tan\phi\sin\left(\psi p + \psi T\right) + \cos\left(\psi p + \psi T\right))}{(\tan\phi\sin\left(\psi p + \psi T\right) + \cos\left(\psi p + \psi T\right))} = \frac{P\mathbf{1}(1 - \tan\phi p \tan\phi - W\mathbf{1}(\tan\phi p \cos\psi p - \sin\psi p))}{(\tan\phi\sin\left(\psi p + \psi T\right) + \cos\left(\psi p + \psi T\right))}$$

The magnitude of *Ts*, is deduced as follows:

$$T\mathfrak{r} = \frac{P\mathfrak{1}(\mathfrak{1} - \tan\phi p \tan\phi - W\mathfrak{1}(\tan\phi p \cos\psi p - \sin\psi p))}{(\tan\phi \sin\left(\psi p + \psi T\right) + \cos\left(\psi p + \psi T\right))}$$

In general, as assumed before that the angle friction for the interfaces of the rock block columns denoted ð Þ *ϕd* is the same as the friction angle on the bases plane ð Þ *ϕ*p, *ϕc* , hence *ϕd* ¼ *ϕp* ¼ *ϕc* ¼ *ϕ*.

Therefore *Ts* is resolved using Eq. (33).

$$\begin{aligned} Ts = \frac{P\mathbf{1}(\mathbf{1} - \tan^2 \phi) - W\mathbf{1}(\tan \phi \cos \psi p - \sin \psi p)}{(\tan \phi \sin (\psi p + \psi T) + \cos (\psi p + \psi T))} \end{aligned} \tag{32}$$

However, if the counter-tilted failure plane dipping at angle *ψc* exists considering the horizontal line, P*Fx* and P*Fy* are solved using Eqs. (33) and (34) respectively.

$$\sum F\mathbf{x} = \mathbf{0} = Pn - Pn + \mathbf{1} - Wn\sin\varphi p + \mathbf{S}n\tag{33}$$

$$\sum Fy = 0 = Rn + Pn\tan\phi - Wn\cos\psi p - Pn + 1\tan\phi\tag{34}$$

If the counter-tilted failure plane at angle *ψc* exists considering the horizontal line, then, when the force T is applied to block 1, the normal force ð Þ *R*1 and shear force ð Þ *S*1 on the base plane of the block 1 are reduced as follows;

$$R\mathbf{1} = P\mathbf{1}\tan\phi + T\sin\left(\wp c + \wp T\right) + W\mathbf{1}\cos\psi p$$

$$\mathbf{S1} = P\mathbf{1} - T\cos\left(\wp c + \wp T\right) + W\mathbf{1}\sin\psi c$$

By resolving the *Fy* equation for *Rn* and making substitutions into the *Fx* equation after counter-tilting of the failure plane along the counter-tilted weak plane *ψc* with consideration of the horizontal line, hence, we set the moments about the pivot point 0 to zero as given below;

$$\begin{aligned} 0 &= P\mathbf{1}(1 - \tan\phi c \tan\phi - W\mathbf{1}(\tan\phi c \cos\psi c - \sin\psi c) - Ts(\tan\phi \sin\psi c + \psi T) \\ &+ \cos\left(\psi c + \psi T\right)Ts(\tan\phi \sin\left(\psi c + \psi T\right) + \cos\left(\psi c + \psi T\right) = P\mathbf{1}(1 - \tan\phi c \tan\phi) \\ &- W\mathbf{1}(\tan\phi c \cos\psi c - \sin\psi c) \end{aligned}$$

We resolved *Ts* as follows:

$$\frac{\operatorname{Ts}(\tan\phi\sin\left(\psi c + \psi T\right) - \cos\left(\psi c + \psi T\right))}{(\tan\phi\sin\left(\psi c + \psi T\right) + \cos\left(\psi c + \psi T\right))} = \frac{P\mathbf{1}(1 - \tan\phi c \tan\phi - W\mathbf{1}(\tan\phi c \cos\psi c - \sin\psi c))}{(\tan\phi\sin\left(\psi c + \psi T\right) + \cos\left(\psi c + \psi T\right))}$$

*Ts*, is deduced as follows:

$$\text{T5} = \frac{P\mathbf{1}(\mathbf{1} - \tan\phi c \tan\phi - W\mathbf{1}(\tan\phi c \cos\psi c - \sin\psi c))}{(\tan\phi \sin\left(\psi c + \psi T\right) + \cos\left(\psi c + \psi T\right))}$$

Based on the earlier assumptions of the angle of friction of the interfaces of the rock columns denoted ð Þ *ϕd* being equal to the angle of friction on the base plane ð Þ *ϕ*p, *ϕc* , thus *ϕd* ¼ *ϕp* ¼ *ϕc* ¼ *ϕ*.then *Ts* is resolved based on the Eq. (35)

$$Ts = \frac{P\mathbf{1}(\mathbf{1} - \tan^2 \phi) - W\mathbf{1}(\tan \phi \cos \psi c - \sin \psi c)}{(\tan \phi \sin \left(\psi c + \psi T\right) + \cos \left(\psi c + \psi T\right))}\tag{35}$$

Stabilization of rock slope by cable forces using limit equilibrium principles demonstrates that support of the "Keystone" is remarkably effective in increasing the scores of stability of the rock slopes prone to toppling. Supporting the rock block around the slope's toe under the influence of toppling failure just at limit equilibrium is simply metastable and its metastability depends on the details of the geometric arrangement of the toppling blocks. On the other hand, reducing the strength of the "keystone" of a slope under toppling that is near failure leads to severe problems.

#### **3. A practical application: Toppling failure**

Failure by toppling of rock columns (**Figure 5**) occurred in shale rock fomations (Shale with grit, Shale with grit type 1 and Shale with grit type 2 & 3) on the northwall of the Nchanga Copper and Cobalt Open Pit situated in the mining city of Chingola in Zambia in June, 2016 [1, 4, 9, 10]. In terms of geology, the Nchanga geology is majorly controlled by Nchanga syncline regional structure with an east–west trend and at 5°-7° plunge to the north. The south limb of the Nchanga Open Pit is noted to be shallow and

**Figure 5.** *Toppling failure of the slope composed of a counter-tilted failure plane.*

#### *Slope Stability Analyses Subjected to Slide Head Toppling Failure Mechanisms DOI: http://dx.doi.org/10.5772/intechopen.108762*

dips to the north at 25°-35° in general together with local changes of shallow to steep dips unlike the north limb which dips at 60°. The Nchanga open pit sub-surface geology is composed of the upper and lower roan dolomites which are usually dipping neaerly horizontal and a series of sedimentary rocks formations overlays them through to basal conglomerates that dips at an angle of 60° towards the north. The slope of the north wall is cut in the ore-body northwall lithology namely; Dolomitic Schist Shale with Grit, Banded Shale, Shale with Pink Quartzite, Grit (type 1,2 & 3), upper and lower roan dolomite Chingola Dolomite, Arkose, Feldspathic Quartzite, Banded Sandstone, Basal Conglomerates partly illustrated in **Figure 5** for the geological engineering profile. Feldspathic Quartzite formation hosts the copper and cobalt that are being mined at Nchanga Open Pits. A summary of the general geological stratigraphy of the mine is shown in **Table 1**. The general slope is elevated at 1330 m – 880 m above mean sea level (asl) from the top and bottom of the elevataions respectively. The overall slope angle is 40° with an overall height of 450 m. An asymmetrical syncline limb on the north wall, the rock columns of the three rock formations (shale with grit type 1, shale with grit type 2, shale with grit) of the rock slope with counter-tilted weak surface overturned. The slope with counter-tilted weak surface under study is sited on the elevations ranging between 1170 m- 970 m above mean sea level (ASL).


#### **Table 1.**

*The general geological layout and rock characteristics.*

Before slope failure, the zone under study being mined at, 155 m high (H), and at a slope angle of 65° (*ψf* ) in a layered shale rock mass (shale with grit type 1, shale with grit type 2shale with grit) dipping at 60° into the face (*ψd* = 60°); with the width (Δ*x*) of individual rock column at 10 m. The angle above the surface of the pit slope was determined to be 5° (*ψs*), and the base of the rock block columns were stepped at b = 1.0 m. Cracks of less than 1 cm on the weak surface started to develop and were noticed when the pit was excavated to a depth of 145 m. The observed weak plane surface on the upper part of the slope (crest) was predicted to daylight on the slope's

**Figure 6.**

*(a) Limiting equilibrium analysis of a block toppling on a rock slope- conceptual slope set up and its geometry. (b) Limiting equilibrium analysis of a block toppling on a slope: Variation/distribution of normal forces (R) and shear (S) forces on the bases of the blocks.*

*Slope Stability Analyses Subjected to Slide Head Toppling Failure Mechanisms DOI: http://dx.doi.org/10.5772/intechopen.108762*

face when the pit was to be mined at 250 m deep since the planned pit depth was at 450 m. No essential remedial measures were therefore undertaken against slope failure and no disruption of the mining operations in the pit was noted. However, based on the geological variations of the Shale rock formations at Nchanga Open Pit, failure surface inclined at 45° on the slope's upper most section day-lighted on the slope's face when the pit was at 144 m deep. It is well appreciated that the plane of weak surface dipping at an angle of 45° in shale with grit rock formation underwent counter-tilting to an angle of 35° in Shale with grit types 1 and 2 formations based on the variations in geological characteristics.

#### **3.1 Calculational procedure**

With regard to the geometry of the slope as given in **Figure 5**, it was presumed that, therefore *ψb* =55°. Using Eqs. (4), n = 14 see **Figure 6a** and rock block column 8 is at the slope's top most part (crest). Using Eqs. (5)–(7), the constants calculated are a1 = 5.0 m, a2 = 5.0 m and b = 1.0 m. These constants are then applied to determine the height *yn* of individual blocks, as well as the height to width ratio. The shale with grit type 1 and shale with grit rock formations' unit weights at Nchanga were noted to be 23kN/m3 3 and 25kN/m respectively. It is assumed that there are no noticeable external forces and that the friction angles are *ϕp* ¼ *ϕd* =30° and *ϕc* =30°. The base plane dips at angles of *ψp* =45° and *ψc* =35°. Therefore, cot *ψp* = 1.0, hence, block columns 14, 13 and 12 are noted to be stable, because for each block indicated, the height to width ratio is noted to be less or equal to1.0. In general, these three blocks are short and their center of gravity was observed to lie inside the base of the block. Block 11 topples because the height-width ratio has a value 1.6, which is greater than 1.0. Therefore, *P*11 is equal to 0 and *P*10 is determined as the greater of *P*10,t and *P*10,*s* given by Eqs. (20)–(23) respectively.


#### **Table 2.**

*Limiting equilibrium analysis of a columnar block toppling slope- listing block calculated forces, stability modes and dimensions.*

This analytical procedure is key for stability determination of individual blocks, based on the failure surface angle of *ψp* moving downwards up to n = 6. At block denoted as n = 5 the original failure plane is counter tilted at angle *θr*, then *ψc* is substituted for *ψp* to determine the stability of individual blocks in turn moving down until n = 1 see Eqs. (20)–(23). The obtained block dimensions, the calculated forces and the stability modes are listed in **Table 2** which indicates that Pn � 1,t is the larger of the two forces up until a value of *n* = 5, where upon Pn � 1,s is larger. Hence, blocks 6 to 11 was noted to be the potential toppling zone, and blocks 1 to 5 denoted a sliding zone. The factor of safety for this slope is determined by varying (increasing) the friction angles until the base plane blocks are stable. It is noted that the required friction angle for limit equilibrium conditions to be satisfied is 36°, or 0.96 (tan 35*/*tan 36). If tan *ϕ* is reduced to 0.577, blocks 1 to 5 in the region around the toe will slide while blocks denoted as 6 to 11 will topple. The anchor tension installed at an angle of 25° through block 1, required to maintain equilibrium, is resolved to be 100kN/m of slope toe based on the Eqs. (24)–(27). This is not a big number, which demonstrates that support of the "Keystone" produces very effective results in increasing stability. On the other hand, reducing the strength of the "keystone" of a slope under toppling, nearing failure, leads to serious consequences. With the definition of the distribution of P forces in the sliding region, the forces Rn and Sn on the base plane of the base blocks is calculated using Eqs. (18) and (19). With the following assumption [*Qn* � 1 ¼ *Pn* � 1 tan *ϕs*], the forces Rn and Sn are determined the region of sliding. **Figure 6b** illustrates the variation of the forces throughout the whole slope. The conditions defined by Rn > 0 and Sn| < *Rn* tan *ϕp* are satisfied everywhere.

#### **4. Summary of the chapter**

A presentation of the mathematical (analytical) method for determination of the counter-tilted failure plane angle for the discontinuous rock slope subjected to block toppling failure mechanism based on the searching technique has been presented. Geometry parameters forming the systematic arrangement of the jointed rock slope that maybe influenced by the presence of the failure planes such as; the application points for the shear and normal forces that acts on the bases of the weak planes and the number of rock columns, the overall base angles of the weak plane angles, forming the regular system have been developed. The modified limit equilibrium method for quantifying the failure mechanisms in the jointed rock slopes having a counter-tilted failure plane prone to block toppling has been proposed. This involves an iterative process in which the dimensions of all the individual blocks and the existing forces acting on them are calculated, and then each block's stability is examined, starting at the uppermost block. Each block will be stable, topple or slide, with the overall stability of the slope being rendered unstable if the lower most block either slides or topples. The iterative process involves determinations of sliding and toppling of the nth block with and without consideration of counter-tilted weak plane respectively. The proposed method of analysis has the potential to curtail errors incurred due to the previous assumptions of the single weak plane for quantifying the failure mechanisms of toppling rock slopes in physical situations with two planar weak planes. It is further noticed that the existence of counter-tilted failure plane in the discontinuous rock slope influence the geometry parameters such as; overall base angles, number of blocks that forms the regular systematic arrangement of the jointed rock slope and application positions *Mn* and *Ln* and the sliding and toppling forces.

*Slope Stability Analyses Subjected to Slide Head Toppling Failure Mechanisms DOI: http://dx.doi.org/10.5772/intechopen.108762*

#### **Author details**

Victor Mwango Bowa<sup>1</sup> \* and Abdul Samson1,2

1 Mining Engineering Department, School of Mines and Mineral Sciences, Copperbelt University, Kitwe, Zambia

2 Malawi University of Business and Applied Sciences, Malawi

\*Address all correspondence to: victor.bowa@cbu.ac.zm

© 2022 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

### **References**

[1] Bowa VM, Xia YY. Influence of counter-tilted failure surface angle on the stability of rock slopes subjected to block toppling failure mechanisms. Bulletin of Engineering Geology and the Environment. 2018a;**77**(4):120-140

[2] Cai JS. Mechanism Research of Toppling Deformation for Homogeneous Equal Thickness Anti-Dip Layered Rock Slopes [Thesis]. Beijing: China University of Geosciences; 2013

[3] Alejano LR. Complex rock slope failure mechanisms involving toppling. In: XIV International Conference about Slope Stability (2018) April, 10–13. Seville, Spain; 2018

[4] Bowa VM, Xia YY. Modified analytical technique for block toppling failure of rock slopes with counter tilted failure surface. Indian Geotechnical Journal. 2018b;**48**(4):713-772

[5] Adhikary DP, Dyskin AV, Jewell RJ, Stewart DP. A study of the mechanism of flexural-toppling failure of rock slopes. Rock Mechanics and Rock Engineering. 1997;**30**:75-93

[6] Pritchard MA, Savigny KW. Numerical modeling of toppling. Canadian Geotechnical Journal. 1990;**27**: 823-834

[7] Pritchard MA, Savigny KW. The Heather Hill landslide: An example of a large scale toppling failure in a natural slope. Canadian Geotechnical Journal. 1991;**28**:410-422

[8] Bowa VM, Wenping Gong WP. Analytical technique for stability analyses of the rock slope subjected to slide head toppling failure mechanisms considering groundwater and

stabilization effects. International Journal of Geo-Engineering. 2021;**2021** (12):4

[9] Bowa VM, Xia Y. Stability analyses of jointed rock slopes with counter-titled failure surface subjected to block toppling failure mechanisms. Arabian Journal for Science, Research Article: Civil Engineering. 2018;**43**(10):5315- 5331

[10] Bowa VM, Xia Y, Kabwe E. Analytical technique for stability analyses of rock slopes subjected to block toppling failure. Mining Technology. 2019;**127**(4):219-229

#### **Chapter 4**

## Numerical Investigation of Rainfall Infiltration-Induced Slope Stability Considering Water-Air Two-Phase Flow

*Wenjing Tian, Herman Peiffer, Benny Malengier, Gang Liu and Qingchao Cheng*

#### **Abstract**

For insights into rainfall infiltration on soil slopes and coupled transmission mechanisms, two-phase flow and finite element analysis were employed to examine water and air movement during the Shuping landslide. The results indicated a division of the landslide surface into two zones: an upper inflow area and a lower overflow area, driven by contrasting inflow and outflow directions. The total water and air flux remained stable, minimally affected by external factors such as rainfall attributes, surface runoff, and air temperature variations. In the inflow area of the slope surface, when rainfall intensity was greater than the total rate of the infiltration of water and air, the magnitude of infiltration equalled to the total rate infiltration of water and air, and runoff generation occurred in this area. Conversely, when infiltration matched rainfall intensity, runoff was absent in this area. In addition, water pressure in the saturated area of the slope surface can be transferred to the groundwater of the slope by pore air pressure, which could also increase the pressure head of the groundwater, and this was also detrimental to slope stability. Regarding uniform rainfall, it significantly reduces the safety factor, potentially making it the most hazardous pattern for slope failure.

**Keywords:** water-air two-phase flow, rainfall infiltration, finite element method, unsaturated seepage, slope stability

#### **1. Introduction**

Rainfall infiltration of a slope is a complex unsaturated seepage process, which is primarily driven by the coupling effect of water and air fluid in the soil [1–3]. The infiltration intensity is influenced by various factors, such as soil permeability, water-air viscosity, and slope boundary conditions [4, 5]. Understanding these factors is essential to predict the rainfall infiltration process and its impact on slope stability [6–8].

Most previous studies [9–12] have used the Richard's rainfall infiltration equation to describe the infiltration process. However, this equation assumes that the pore air pressure is constant and equal everywhere (usually set to atmospheric pressure), which essentially ignores the coupling and cooperation between water and air in the soil. This simplification may be valid under certain circumstances, but it may not fully capture the complex interactions between water and air in real-world engineering scenarios. Numerous research findings [13–17] indicated that, in the case of a closed boundary condition, the effect of pore air on the flow of the liquid phase (water) cannot be neglected. During the rainfall infiltration process, the liquid phase (water) exerts pressure on the air in the soil, resulting in a jacking force that opposes the direction of water movement. This phenomenon can cause a decrease in the rainfall infiltration rate because of the reduction of the water pressure gradient. This effect of air on the infiltration process has been observed and reported in numerous studies [18, 19]. Therefore, the consideration of water-air coupling effect is important for accurate modelling and prediction of rainfall infiltration in soil.

Some researchers have demonstrated that the water and air infiltration process of slope rainfall can be accurately captured by considering the influence of the air phase. Based on their findings, it is recommended to adopt the two-phase flow theory of water and air for a more practical exploration of the problem of slope rainfall infiltration. This chapter employs the basic theory and methodology of water- air two-phase flow to simulate and analyse the rainfall infiltration process of Shuping landslide in the Three Gorges Reservoir area, with the aim of investigating the general law of slope rainfall infiltration. While most studies on two-phase flow primarily focus on exploring changes in pressure and other parameters during infiltration and verifying the impact of the air phase, there is a lack of research on the overall induction of rainfall infiltration laws. By studying the general law of slope rainfall infiltration, this research serves as a complement to existing laws of rainfall infiltration and provides a more comprehensive basis for engineering applications. In addition, the effect of rainfall pattern on slope failure is also analysed in this research.

#### **2. Governing differential equation for water-air two-phase flow**

Liquid phase (water) flow equation and air phase flow equation are included in the governing differential equations for water-air two-phase flow. These two equations can realize the coupling process by correlation of some variables including saturation, pore air pressure, matric suction as well as porosity. The flow of water and air not only conforms with mass conservation equation driven by pressure but also considers some variation factors including soil deformation, compressibility, viscosity as well as saturation. The governing differential equations for water-air two-phase flow can be described as follows [20]:

$$\mathbf{nS}\frac{\partial \ln \rho\_l}{\partial t} + \mathbf{S}\frac{\partial n}{\partial t} + n\frac{\partial \mathbf{S}}{\partial t} + \nabla \cdot \left[ -\frac{k\_r^l k}{\mu\_l} \cdot \left( \nabla p\_l + \rho\_l \mathbf{g} \right) \right] - \frac{\mathbf{1}}{\rho\_l} \mathbf{Q}\_l = \mathbf{0} \tag{1}$$

$$\ln(\mathbf{1} - \mathbf{S}) \frac{\partial \ln \rho\_{\mathbf{g}}}{\partial t} + (\mathbf{1} - \mathbf{S}) \frac{\partial \mathbf{n}}{\partial t} - n \frac{\partial \mathbf{S}}{\partial t} + \nabla \cdot \left[ - \frac{k\_r^{\mathbf{g}} \mathbf{k}}{\mu\_{\mathbf{g}}} \cdot \left( \nabla p\_{\mathbf{g}} + \rho\_{\mathbf{g}} \mathbf{g} \right) \right] - \frac{\mathbf{1}}{\rho\_{\mathbf{g}}} \mathbf{Q}\_{\mathbf{g}} = \mathbf{0} \tag{2}$$

*Numerical Investigation of Rainfall Infiltration-Induced Slope Stability Considering… DOI: http://dx.doi.org/10.5772/intechopen.113723*

Where n represents the porosity; *Sr* represents the degree of saturation; *ρβ* represents the density of the phase (β ¼ 1 represents the water phase, β ¼ grepresents the air phase); *Q<sup>β</sup>* represents the source of phase (kg/m); k represents the intrinsic permeability of soil determined by pore characteristics (*m*2); *k<sup>β</sup> <sup>r</sup>* represents the relative permeability coefficient of the phase; *μβ* represents the viscosity coefficient of the phase N � <sup>s</sup>*=*m<sup>2</sup> ð Þ; *<sup>p</sup><sup>β</sup>* is the pore pressure of the <sup>β</sup> phase N*=*m<sup>2</sup> ð Þ; g is the gravitational acceleration (N/kg).

The dynamics of water-air two-phase flow are described by a complex system of partial differential equations that are nonlinear and depend on both space and time variables. To tackle the spatial domain, the Galerkin finite element method is utilized, while the time domain is discretized using a difference method. This approach enables a robust and accurate calculation of the flow dynamics by considering two unknowns, namely the pore pressure (ρ*g*) and the saturation (S*r*), which are solved iteratively through a cyclic calculation of Eq. (1) and Eq. (2). This iterative process ensures that the solution converges to the correct values. Overall, this methodology provides an efficient and reliable approach for investigating water-air two-phase flow dynamics, which can help to understand the behaviour of fluid flow in porous media, and provide insights into the optimization of industrial processes.

#### **3. Computational model**

In pursuit of a more profound understanding of the underlying principles governing rainfall infiltration, this section employs the illustrative framework of the Shuping landslide case study. Through the application of advanced simulation techniques, the dynamic interplay between water and air within the slope during the process of rainfall infiltration is meticulously recreated. This simulation endeavour is strategically designed to shed light on the fundamental tenets that dictate the behaviour of rainfall infiltration on slopes. It takes into account the distinctive attributes of this process, the variations in its intensity, and its propensity for generating runoff.

Subsequent to this, the discussion takes a comprehensive turn by delving into the intricate mechanisms that underpin the intermingling of water and air within the environmental context of the slope. By deciphering these mechanisms, the study aims to offer deeper insights into the multifaceted factors influencing the behaviour of water and air during rainfall infiltration. The discourse seeks to unravel the dynamics that govern their coexistence within the slope's environment, ultimately contributing to a more holistic and accurate comprehension of this intricate natural process.

#### **3.1 Study area and geometric model**

This study presents a case study of the Shuping landslide, which is located on the north bank of the Yangtze River, 61 km away from the dam site of the Three Gorges Project. The landslide body comprises a loose rock clastic layer and a loess soil layer, which have good permeability and are conducive to groundwater infiltration. The bottom of the landslide body is composed of phyllite, which has relatively weak permeability and can form a relative waterproof layer.

The Shuping landslide is characterized as a relict landslide, characterized by a historical chronicle of recurring sliding occurrences. Since the initiation of the Three Gorges Reservoir impoundment in June 2003, a persistent sequence of deformations

**Figure 1.** *Finite element grid diagram of the slope cross-section.*

has been discerned. Consequently, there has been a progressively escalating manifestation of surface deformations within the sliding mass. Particularly notable is the observable trend of continuous propagation observed in the emergence of cracks within the slope. Furthermore, instances of localized significant collapses have been identified. Presently, the landslide remains subjected to an ongoing process of deformation.

The Shuping landslide boasts a historical record marked by recurrent events of collapse and sliding. Its stability is relatively compromised, thereby rendering it susceptible to deformations and failures instigated by external variables such as precipitation and fluctuations in reservoir water levels. If a substantial decrease in the Three Gorges Reservoir water level, reducing it from 175 meters to 145 meters, coincides with intense rainfall, the sliding mass would be susceptible to incurring substantial deformations, thereby possibly precipitating a precarious state of instability.

The thickness of the sliding body provides the material foundation for the landslide. The front edge of the Shuping landslide is 1100 m wide, the back edge is 800 m wide, and the longitudinal length is 550 m. The back edge is distributed in the line of the elevation of 280–338 m. For numerical modelling, the main sliding section is considered, with a horizontal length of 600 m and a vertical height of 350 m. To accurately simulate rainfall infiltration, the surface layer of the landslide body is discretized into elements with varying thicknesses, ranging from 0.25 to 0.75 m along the slope depth direction. The finite element calculation grid obtained by subdivision consists of 3499 elements and 3427 nodes, as shown in **Figure 1** (geometric model for numerical simulation). This approach provides a comprehensive understanding of the landslide characteristics and the factors that influence its behaviour, which can contribute to the development of effective mitigation strategies for similar geological hazards.

#### **3.2 Initial and boundary conditions**

To reduce the uncertainty impact of initial state on this research, this study incorporates measured rainfall data from the Three Gorges Reservoir area spanning several years (provided by the Headquarters of Geological Disaster Prevention in the Three Gorges Reservoir Area of the Ministry of Land and Resources). The seepage

#### *Numerical Investigation of Rainfall Infiltration-Induced Slope Stability Considering… DOI: http://dx.doi.org/10.5772/intechopen.113723*

field of the slope is calculated and simulated in the fully saturated state until it reaches to a relatively stable state, which is used as the initial condition for subsequent calculations. By adopting this approach, the study aims to minimize potential errors or biases that may arise from using arbitrary initial conditions. The use of measured data and stable conditions provides a more accurate and realistic representation of the initial state of the system, which is critical for achieving meaningful and reliable results. This study conducts a thorough examination of the real-world system under investigation. It investigates various aspects, including the behaviour of materials, the system's boundaries, its geometries, and significant variables such as the impact of rainfall. The interplay of these diverse elements collectively influences the system's behaviour. Through a detailed analysis of these factors, the study aims to unravel the intricate dynamics that govern the system's behaviour within its environmental context.

In this research, the side and bottom of the rear edge of the model are treated as impervious boundaries, and the water level and its variation at the boundary are calculated. The slope surface and leading edge below the reservoir water level are considered as known water pressure boundaries, and their values are dependent on the reservoir water level elevation (145 m). On the other hand, the slope surface above the reservoir water level is treated as the known air pressure boundary. As the runoff water head on the slope is relatively small compared with the atmospheric pressure, it is assumed to be negligible. In the present analysis, the air pressure on the air pressure boundary is set to atmospheric pressure. The proposed approach enables the determination of infiltration rates of both water and air from calculation results, without the need for their *a priori* specification as flow boundaries. Compared to the traditional single-phase flow calculations, this approach provides a more comprehensive understanding of the dynamics of the water-air two-phase flow, which is important for achieving more accurate predictions and more effective risk assessments in geotechnical engineering applications.

#### **3.3 Constitutive relations and parameters**

For the given slope, the mathematical description of the rainfall infiltration process is determined by the governing differential equation, initial and boundary conditions, constitutive relation, and corresponding parameters of water and air two-phase flow. Eq. (1) and Eq. (2) involve five unknowns: S, *k<sup>l</sup> <sup>r</sup>*, *<sup>k</sup><sup>g</sup> <sup>r</sup>*, *pl* and *pg*, which are solved concurrently by incorporating the relationship between soil-water characteristics determined by soil properties and the relative permeability function of water and air.

In this research, the commonly used Van Genuchten model is adopted to represent the soil-water characteristic curve. The mathematical relationship between matric suction *pc* � *pg* � *pl* and saturation is expressed as follows [21]:

$$p\_{\varepsilon} = -p\_0 \left( \mathbf{S}\_{\varepsilon}^{-1/\lambda} - \mathbf{1} \right)^{1-\lambda} \tag{3}$$

Where S*<sup>e</sup>* represents the effective water saturation *Se* <sup>¼</sup> *<sup>S</sup>*�*Srl* <sup>1</sup>�*Srl* ; *Srl* represents the residual water saturation; *p*<sup>0</sup> and λ are parameters of this model. According to some literature data [22], this calculation adopts the value of *p*<sup>0</sup> ¼ 1*:*33, *λ* ¼ 0*:*41, *Srl* ¼ 0*:*15 as three parameters.

The relative permeability coefficients of water (*k<sup>l</sup> r*) and air (*kg <sup>r</sup>*) are using the Van Genuchten-Mualem [22] model and Corey model, respectively. In this study, the relative permeability coefficients that contain effective saturation variable can be expressed as follows [23, 24]:

$$k\_r^l = \sqrt{S\_\epsilon} \left\{ \mathbf{1} - \left[ \mathbf{1} - (\mathbf{S}\_\epsilon)^{1/\lambda} \right]^\lambda \right\}^2 \tag{4}$$

$$k\_r^\mathbf{g} = (\mathbf{1} - \mathbf{S}\_\epsilon)^2 \left[\mathbf{1} - \mathbf{S}\_\epsilon^2\right] \tag{5}$$

The values of some other parameters are as follows [22]: n ¼ 0*:*15, <sup>k</sup> <sup>¼</sup> <sup>4</sup>*:*<sup>0</sup> <sup>∗</sup> <sup>10</sup>�12*m*2, g <sup>¼</sup> <sup>9</sup>*:*8N*=*kg, *<sup>ρ</sup><sup>g</sup>* <sup>¼</sup> <sup>1</sup>*:*29*kg=m*3, *<sup>ρ</sup><sup>l</sup>* <sup>¼</sup> <sup>1</sup> <sup>∗</sup> <sup>10</sup><sup>3</sup> *kg=m*3, *<sup>μ</sup><sup>l</sup>* <sup>¼</sup> <sup>1</sup> <sup>∗</sup> <sup>10</sup><sup>3</sup> *<sup>N</sup>* � *<sup>s</sup>=m*2, *<sup>μ</sup><sup>g</sup>* <sup>¼</sup> <sup>1</sup> <sup>∗</sup> <sup>10</sup>�<sup>5</sup> *<sup>N</sup>* � *<sup>s</sup>=m*2.

#### **3.4 Slope stability analysis method considering pore air pressure**

In this study, the slope stability analysis is conducted to use the residual thrust method, which incorporates the consideration of matric suction and pore air pressure. The calculation procedure involves several steps: an initial assumption of the safety factor, subsequent calculation of the thrust starting from the first slide at the top of the slope and progressing towards the last slide, and finally determining the safety factor value at which the thrust reached zero. This zero-thrust condition represents the equilibrium state, thus yielding the final safety factor value [25–28].

$$\begin{aligned} P\_i &= W\_i \sin \alpha\_i - \frac{c\_i' l\_i + \left(W\_i \cos \alpha\_i - p\_a' l\_i\right) \tan \phi\_i' + l\_i \left(p\_a^i - p\_w^i\right) \tan \phi^b}{F\_s} + P\_{i-1} \psi\_i \\ \psi\_i &= \left[\cos(a\_{i-1} - a\_i) - \frac{\tan \phi\_i'}{F\_s} \sin(a\_{i-1} - a\_i)\right] \end{aligned} \tag{6}$$

Where *Fs* is the safety factor of the slope, *Pi* is the sliding force of soil slice, *c*<sup>0</sup> *<sup>i</sup>* is the effective cohesive force of slice, *φ*<sup>0</sup> *<sup>i</sup>* is the effective internal friction angle of slice, *li* is the width of soil slice, *Wi* is the weight of soil slice, *α<sup>i</sup>* is the angle of bottom soil slice, *pi <sup>a</sup>* and *pi <sup>w</sup>* are pore air pressure and pore water pressure of slice, *tan φ<sup>b</sup>* is the rate of shear strength that increases with an increase in matric suction, and *ψ<sup>i</sup>* is the transfer coefficient of slice i. The model recommended by Vanapalli [29] in this research can be expressed as a function of effective saturation. The equation is as follows:

$$
\tan \, q^b = \mathbb{S}^\* \cdot \tan \, q' \tag{7}
$$

Where the effective saturation is expressed as follows:

$$\mathcal{S}^\* = (\mathcal{S}\_r - \mathcal{S}\_{rw})/(\mathbf{1} - \mathcal{S}\_{rw}) \tag{8}$$

#### **4. Rainfall infiltration analysis of slope**

In this research, we present the results of finite element calculations of two-phase flow (water and air) in Shuping landslide. Our findings reveal the general rule of rainfall infiltration in slopes, examining characteristics such as infiltration intensity,

*Numerical Investigation of Rainfall Infiltration-Induced Slope Stability Considering… DOI: http://dx.doi.org/10.5772/intechopen.113723*

flow-producing strips, and infiltration patterns. Additionally, we elucidate the coupling force transfer mechanism between water and air in slopes. Through our analysis, we aim to contribute to a deeper understanding of the complex processes that govern the behaviour of landslides under the influence of rainfall infiltration.

#### **4.1 Rainfall infiltration characteristics of slope surfaces**

According to the calculation results of Shuping landslide, the direction and total velocity of material (water, air) entering and leaving the slope surface show a certain regularity. According to different directions of material (water and air) inflow and outflow, the slope surface can be divided into inflow area and overflow area. If the fluid tends to flow in the direction of the slope, exhibiting an upward or uphill movement, this area is commonly known as the inflow area. Conversely, the overflow area represents the zone where the flow rate of water or air predominantly moves away from the slope (**Figure 2**). The distribution of inflow area and overflow area of slope surface has strong regularity. The water and air inflow area is usually located in the middle and upper part of the slope, while the overflow area is generally located in the lower part. The boundary line between the suction zone and the overflow zone is usually not fixed. It is not only related to the geometric characteristics and permeability characteristics of slope, but also affected by rainfall characteristics such as rainfall duration, rainfall interval, and rainfall intensity. With the extension of rainfall duration, the boundary line between inflow area and overflow area will move upward, the inflow area of slope surface will shrink, and the overflow area will increase. With the increase of rainfall interval, the boundary line will gradually move down, and the inflow area of slope surface will increase, while the overflow area will decrease.

Based on numerical analysis, the total flow rate of water and air through the slope surface exhibits phase-pair stability, and its magnitude is primarily influenced by slope geometry and the transmission characteristics of geotechnical materials, while external factors such as rainfall characteristics, slope flow production, and air temperature variation have relatively little effect. The flow rate of water and air into and out of the slope surface provides a measure of the intensity of rainfall inflow and overflow. **Figure 3** depicts the relationship between the intensity of water and air inflow (overflow) and time. The curve indicates that, as the rainfall duration increases, the maximum inflow and overflow intensities of the Shuping landslide slope attain a certain stability, with the

**Figure 2.** *Slope surface (water, air): Inflow and overflow area.*

**Figure 3.** *Variation of inflow and overflow intensity over time.*

former being around 0.4 mm/min and the latter around 0.5 mm/min. The average inflow and overflow intensity values of water and air on the slope exhibit a slight downward trend, with both values being quite similar at approximately 0.1 mm/min.

#### **4.2 Rainfall infiltration intensity and flow-producing conditions on slope surface**

The assessment of rainfall infiltration intensity and runoff production on slope surfaces is closely linked to the relative magnitude of infiltration and the total inflow rate of the slope surface, comprising both water and air. In the inflow area of the slope surface, infiltration of water and air relies primarily on the composition of water and air at the slope's boundary. In the absence of rainfall on the surface, air is drawn into the inflow area through its boundary. Conversely, when rainfall occurs on the surface, if the rainfall intensity is lower than the total inflow rate of the slope (consisting of water and air), the slope's inflow area will not produce flow, and the infiltration intensity of rainwater will match the rainfall intensity. The slope's boundary will absorb a mixture of water and air with varying compositions. However, if the rainfall intensity exceeds the total inflow rate of the slope (comprising both water and air), runoff will occur on the slope, and the infiltration intensity will equal the total inflow intensity of water and air. In such cases, the surface of the slope is covered by air, while water is drawn at the inflow area's boundary.

Regarding the Shuping landslide, rainfall intensity is much greater than the total inflow intensity of the slope, approximately 0.1 mm/min, which is a prerequisite for runoff to be generated on the slope. The maximum inflow area is located at the top of the slope, when the rainfall intensity exceeds the maximum inflow intensity of the slope surface, which is approximately 0.4 mm/min, and the slope may experience full flow generation. These observations suggest that rainfall intensity and the total inflow rate of the slope surface are essential factors to consider when analysing the potential for runoff generation and slope stability during rainfall events.

In the overflow area of the slope surface, the overflow of water and air is closely linked to the water and air content at the slope's boundary and the slope saturation. According to a report by the Yangtze River Scientific Research Institute in 2017, when rainfall infiltration intensity was 0, the slope surface produced runoff, and the outflow at the boundary of the overflow area was water. In an unsaturated state, the infiltration intensity of rainfall is equivalent to the intensity of the rainfall itself, resulting in a lack

*Numerical Investigation of Rainfall Infiltration-Induced Slope Stability Considering… DOI: http://dx.doi.org/10.5772/intechopen.113723*

of flow on the slope surface. The outflow from the boundary of the overflow zone consists of a combination of water and air. In scenarios where the entire overflow area is dry, the outflow at the boundary of the overflow area exclusively consists of pore air.

These observations yield insights into the pivotal roles played by several key factors in the process of water and air overflow on the slope surface. Among these factors, the water and air content existing at the boundary of the overflow region, as well as the saturation state of the slope, emerges as critical determinants. Moreover, these findings indicate that a confluence of variables, namely rainfall intensity, slope saturation, and the water and air content at the slope's boundary, collectively governs the production of runoff and the characteristics of overflow events.

In essence, these insights underscore the intricate interplay of factors that orchestrate the dynamics of water and air behaviour within the context of slope surface overflow. The water and air content at the boundary of the overflow zone, alongside the slope's saturation state, serves as foundational considerations that shape the overflow process. Furthermore, the interrelation between these variables and external factors such as rainfall intensity adds layers of complexity to the observed behaviours. These findings contribute to an enhanced understanding of the multifaceted mechanisms governing the overflow dynamics, thereby facilitating more accurate assessments and predictions in geotechnical and hydrological investigations.

#### **4.3 Water-air coupling force transfer mechanism**

Rainfall infiltration is a process where water and air interact within soil pores, influencing each other's movement. Water moves downward due to gravity, causing air to move within the pores. The compression of air generates a counteracting force against the water's movement. However, due to its low viscosity (about 1% that of

**Figure 4.** *Mechanism of water-air coupling in force transfer.*

water), air's effect on slope stability is unfavourable. **Figure 4** depicts the force transfer mechanism resulting from water-air coupling. This coupling leads to pressure on the subsurface water, highlighting the detrimental influence of air on slope stability.

**Figures 5**–**7** show the saturation, pore air pressure, and pore water pressure distribution cloudy map of Shuping landslide at a specific time. It reveals that the saturated water pressure at the top of the slope is transmitted to the saturated area at the slope's front through closed pore air in the unsaturated area. This phenomenon leads to a partial elevation of groundwater level in the anti-sliding area at the slope's front, resulting in the formation of external normal thrust that adversely affects slope stability.

**Figure 5.** *Distribution of water saturation for Shuping landslide at one time.*

#### **Figure 6.**

*Distribution of pore air pressure for Shuping landslide at one time.*

**Figure 7.** *Distribution of pore water pressure for Shuping landslide at one time.*

*Numerical Investigation of Rainfall Infiltration-Induced Slope Stability Considering… DOI: http://dx.doi.org/10.5772/intechopen.113723*

#### **4.4 Effect of rainfall type on infiltration**

By analysing pertinent data, the rainfall pattern can be categorized into four distinct conditions, as depicted in **Figure 8**. These conditions are characterized as four rainfall patterns. Uniform rainfall pattern: during this condition, the rainfall intensity remains relatively constant over the specified time period. Pre-peak rainfall pattern: this condition exhibits a linear decrease in rainfall intensity over time, leading up to a peak point. Post-peak rainfall pattern: in this condition, the rainfall intensity demonstrates a linear increase over time, following a peak point. Mid-peak rainfall pattern: this condition displays a rainfall intensity that initially increases and then subsequently decreases, forming a peak point.

**Figures 9** and **10** illustrate the variations in pore air pressure and safety factor under four different rainfall patterns. In the case of uniform rainfall pattern, when the rainfall intensity remains constant, the saturation of the soil on the slope surface becomes saturated. Consequently, the slope surface forms a relatively impermeable state, which limits the outflow of air from the slope surface. As a result, the pore air pressure continues to rise. With increasing saturation, the matric suction decreases, leading to an increase in pore water pressure. During this period, the rainfall not only increases the load on the slope but also weakens the strength parameters, causing a decrease in the safety factor of the reservoir bank slope. As the rainfall duration increases, the infiltration intensity on the slope surface gradually decreases, resulting in a reduction in the amount of infiltrated water. Consequently, the pore air pressure, pore water pressure, and safety factor eventually reach a stable state. In the case of prepeak rainfall, the initial stage shows a similar pattern to uniform rainfall, where the

**Figure 8.** *Four types of rainfall patterns.*

**Figure 9.** *Temporal variation in pore air pressure associated with four rainfall patterns.*

**Figure 10.** *Temporal variation in safety associated with four rainfall patterns.*

pore air pressure and pore water pressure increase non-linearly, leading to a decrease in the safety factor. During linearly decreasing rainfall, the slope surface gradually transitions from a saturated state to an unsaturated state. Pores in unsaturated soil provide a pathway for the overflow of pore air pressure from the slope body. In comparison with the early stages of rainfall, the pore air pressure decreases rapidly, while the pore water pressure increases in the latter half of the pre-peak rainfall. Consequently, the slope safety factor is gradually recovered. As the rainfall intensity is relatively low in the early period, the pore air pressure begins to increase from 0. As the rainfall continues, the intensity gradually increases. As a result, the surface soil becomes saturated rapidly during the later stage of rainfall, leading to a rapid increase in pore air pressure and pore water pressure, while the safety factor of the slope decreases. In the case of mid- peak rainfall, the rainfall intensity gradually increases during the early stages, and the pore air pressure increases gradually from 0. As a result, the safety factor experiences a slow decrease in the initial stage. With the continuous increase in rainfall intensity, the pore air pressure decreases rapidly, while the pore water pressure increases, leading to a rapid decrease in the safety factor. At the point when the rainfall intensity reaches its maximum, both the pore air pressure and pore water pressure reach their maximum values, resulting in the minimum safety factor. Subsequently, as the rainfall intensity starts to decrease, the pore air pressure and pore water pressure also decrease, causing a slow increase in the safety factor.

*Numerical Investigation of Rainfall Infiltration-Induced Slope Stability Considering… DOI: http://dx.doi.org/10.5772/intechopen.113723*

#### **5. Conclusions**

This study initially presents the fundamental patterns governing rainfall infiltration on slopes. This presentation encompasses various aspects, including slope surface characteristics related to rainfall infiltration, the intensity of rainfall infiltration, and the prevailing runoff conditions. Subsequently, a comprehensive explanation is provided regarding the intricate mechanism underlying the coupled transmission forces of water and air within the slope's environment, which also constitutes the novelty of this work. Furthermore, this research also encompasses an analysis of the impact of rainfall patterns on slope failure.

Based on the principles of water-air two-phase flow, the rainfall infiltration process of Shuping landslide is analysed through finite element method. The following conclusions could be drawn from the analysis:


#### **Acknowledgements**

The work was supported by the Geotechnical Laboratory at Ghent University and Geosound.be. The first author would like to acknowledge the financial support received from the China Scholarship Council (No.201908420298).

### **Conflict of interest**

The authors declare no conflict of interest.

### **Author details**

Wenjing Tian<sup>1</sup> \*, Herman Peiffer<sup>1</sup> , Benny Malengier<sup>2</sup> , Gang Liu<sup>3</sup> and Qingchao Cheng<sup>4</sup>

1 Department of Civil Engineering, Ghent University, Ghent, Belgium

2 Department of Textiles and Chemical Engineering, Ghent University, Ghent, Belgium

3 Department of Hydraulic and Environment Engineering, China Three Gorges University, Yichang, China

4 Hangzhou Guodian Dam Safety Engineering, Hangzhou, China

\*Address all correspondence to: wenjing.tian@ugent.be

© 2024 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

*Numerical Investigation of Rainfall Infiltration-Induced Slope Stability Considering… DOI: http://dx.doi.org/10.5772/intechopen.113723*

#### **References**

[1] Li B, Tian B, Tong FG, Liu C, Xu XL. Effect of the water-air coupling on the stability of rainfall-induced landslides using a coupled infiltration and hydromechanical model. Geofluids. 2022;**8**:1-16. DOI: 10.1155/2022/ 3036905

[2] Liu G, Tong FG, Zhao YT. A force transfer mechanism for triggering landslides during rainfall infiltration. Journal of Mountain Science. 2018;**15**: 2480-2491. DOI: 10.1007/s11629-018- 5043-x

[3] Liu G, Tong FG, Tian B. A finite element model for simulating surface run-off and unsaturated seepage flow in the shallow subsurface. Hydrological Processes. 2019;**33**(26):3378-3390. DOI: 10.1002/hyp.13564

[4] Hou XK, Vanapalli SK, Li TL. Water infiltration characteristics in loess associated with irrigation activities and its influence on the slope stability in Heifangtai loess highland, China. Engineering Geology. 2018;**234**:27-37. DOI: 10.1016/j.enggeo.2017.12.020

[5] Airmo JM, Rahardjo H, Leong EC. Infiltration effects on stability of a residual soil slope. Computers and Geotechnics. 2000;**26**(2):145-165. DOI: 10.1016/S0266352X (99)00035-X

[6] Liu G, Tong FG, Tian B, Tian WJ. Influence of atmospheric temperature on shallow slope stability. Environmental Earth Sciences. 2019;**78**(22):632. DOI: 10.1007/s12665-019-8649-6

[7] Tong FG, Jing L, Zimmerman RW. A fully coupled thermo-hydro-mechanical model for simulating multiphase flow, deformation and heat transfer in buffer material and rock masses. International Journal of Rock Mechanics and Mining

Sciences. 2010;**47**(2):205-217. DOI: 10.1016/j.ijrmms.2009.11.002

[8] Zhou YF, Tham LG, Yan WM, Dai FC, Xu L. Laboratory study on soil behavior in loess slope subjected to infiltration. Engineering Geology. 2014; **183**:31-38. DOI: 10.1016/j.enggeo. 2014.09.010

[9] Koyama A, Fujimoto T, Suetsugu T, Fukubayashi Y. Analysis of embankment failure mechanism in reservoirs due to rainfall infiltration during heavy rainfall. Natural Hazards and Risk. 2022;**13**(1): 1849-1866. DOI: 10.1080/19475705. 2022.2102440

[10] Komori D, Rangsiwanichpong P, Inoue N, Ono K, Watanabe S, Kazama S. Distributed probability of slope failure in Thailand under climate change. Climate Risk Management. 2018;**20**:126-137. DOI: 10.1016/j.crm.2018.03.002

[11] Sharma RH, Nakagawa H. Numerical model and flume experiments of single and two-layered hillslope flow related to slope failure. Landslides. 2010;**7**:425-432. DOI: 10.1007/s10346-010-0205-0

[12] Liu C, Tong F, Li B, et al. A water retention curve model describing the effect of temperature. European Journal of Soil Science. 2020;**71**(1):44-54. DOI: 10.1111/ejss.12825

[13] Sun DM, Zang YG, Semprich S. Effects of airflow induced by rainfall infiltration on unsaturated soil slope stability. Transport in Porous Media. 2015;**107**:821-841. DOI: 10.1007/ s11242-015-0469-x

[14] Sung EC. Stability analysis of unsaturated soil slopes considering water-air flow caused by rainfall

infiltration. Engineering Geology. 2016; **211**:184-197. DOI: 10.101 6/j. enggeo.2016.07.008

[15] Sun DM, Li XM, Feng P, Zang YG. Stability analysis of unsaturated soil slope during rainfall infiltration using coupled liquid-air-solid three-phase model. Water Science and Engineering. 2016;**9**(3):183-194. DOI: 10.1016/j.wse. 2016.06.008

[16] Hu R, Chen YF, Zhou CB. Modelling of coupled deformation, water flow and air transport in soil slopes subjected to rain infiltration. Science China Technological Sciences. 2011;**54**: 2561-2575. DOI: 10.1007/s11431-011- 4504-z

[17] Touma J, Vauclin M. Experimental and numerical analysis of two-phase infiltration in a partially saturated soil. Transport Porous Media. 1986;**1**:27-55. DOI: 10.1007/BF01036524

[18] Xue S, Yang Z, Hu R, et al. Splitting dynamics of liquid slugs at a T-junction. Water Resources Research. 2020;**56**(8): e2020WR027730. DOI: 10.1002/ essoar.10503637.1

[19] Chen P, Wei CF. Numerical procedure for simulating the two-phase flow in unsaturated soils with hydraulic hysteresis. International Journal of Geomechanics. 2016;**16**(1):1532-3641. DOI: 10.1061/(ASCE)GM.1943- 5622.0000505

[20] Tong FG. Numerical Modelling of Coupled Thermo-Hydro-Mechanical Processes in Geological Porous Media. Stockholm, Sweden: KTH; 2010

[21] Van Genuchten MT. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal.

1980;**44**:892-898. DOI: 10.2136/ sssaj1980.03615995004400050002x

[22] Mualem Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resources Research. 1976;**12**:513-522. DOI: 10.1029/WR0 12i003p00513

[23] Corey AT. The interrelation between air and oil relative permeabilities. Producers Monthly. 1954;**19**:38-41

[24] Tian W, Peiffer H, Malengier B, Xue S, Chen Z. Slope stability analysis method of unsaturated soil slopes considering pore air pressure caused by rainfall infiltration. Applied Science. 2022;**12**:11060. DOI: 10.3390/ app122111060

[25] Wei WB, Cheng YM, Li L. Threedimensional slope failure analysis by the strength reduction and limit equilibrium methods. Computers and Geotechnics. 2009;**36**(1-2):70-80. DOI: 10.1016/ j.compgeo.2008.03.003

[26] Fang N, Ji CS, Crusoe GE. Stability analysis of the sliding process of the west slope in Buzhaoba open-pit mine. International Journal of Mining Science and Technology. 2016;**26**(5):869-875. DOI: 10.1016/j.ijmst.2016.05.033

[27] Su AJ, Zou ZX, Lu ZC, Wang JG. The inclination of the interslice resultant force in the limit equilibrium slope stability analysis. Engineering Geology. 2018;**240**:140-148. DOI: 10.1016/ j.enggeo.2018.04.016

[28] Yao W, Li C, Zhan H. Timedependent slope stability during intense rainfall with stratified soil water content. Bulletin of Engineering Geology and the Environment. 2019;**78**: 4805-4819. DOI: 10.1007/s10064-018- 01437-3

*Numerical Investigation of Rainfall Infiltration-Induced Slope Stability Considering… DOI: http://dx.doi.org/10.5772/intechopen.113723*

[29] Vanapalli SK, Sillers WS, Fredlund MD. The meaning and relevance of residual state to unsaturated soils. In: 51st Canadian Geotechnical Conference. Vol 8. Edmonton Alberta, Canada; 1998

#### **Chapter 5**

## The Impact of Earthquakes on Dropout Doline (Cover Collapse Sinkhole) Development: A Case Study from the Environs of Mečenčani and Borojevići (Croatia)

*Márton Veress, Natalija Matić, Zoltán Mitre and Gábor Szunyogh*

#### **Abstract**

In this study, the effect of earthquakes, beginning on 28 December 2020, on dropout doline development in the environs of Mečenčani and Borojevići was investigated. For that purpose, the shape of the doline, the inclination of the bearing surface and the rate of doline development were determined. A further analysis involved the characteristics of groundwater in the environs of the dolines and a functional relationship between the depth and the diameter of the dolines was sought. A model is proposed for the failure of the ceiling of cover cavities without support. The intensity of doline development is explained by favourable environment (dual cavity system, low inclination of the bearing surface, the presence and fluctuation of groundwater, etc.), the direct effect of earthquakes (material failure induced by earthquakes) and by their indirect effect (the partial solifluction of the ceiling material, lowered groundwater level).

**Keywords:** dropout doline, earthquake, groundwater, karst water, cavity

#### **1. Introduction**

In this study, the development of dropout dolines (cover collapse sinkholes) in the environs of Mečenčani and Borojevići (Croatia) to the effect of the earthquakes of Petrinja between 28 December 2020 and 03 March 2021 is interpreted. Croatian researchers compiled a comprehensive documentation on the effect of the Petrinja earthquakes and also described the dolines that were formed during this time [1, 2]. All other investigations are of geomorphological nature [3], those concerned with remediation measures and priorities for immediate action [4], ground displacement using data from orbits of the Sentinel-1 mission action [5, 6]; fault geometry and the coseismic slip distribution [7], earthquake and deformations [8, 9]. To date, not a single investigation

considered the effect of earthquakes on dropout doline development in the environs of Mečenčani and Borojevići such as the shape of the doline, the inclination of the bearing surface and the rate of doline formation. It is very important to emphasize that doline development damaged houses in the villages. The novelty of this investigation includes a model for the failure of the ceiling of cover cavities without support, analysis of the characteristics of groundwater in the environs of the dolines and a functional relationship between the depth and the diameter of the dolines.

Subsidence dolines (subsidence sinkholes) develop on unconsolidated, permeable or partly permeable rock (on concealed karst). Their varieties are dropout dolines (cover collapse sinkholes), suffosion dolines and compaction dolines [10, 11]. Dropout dolines are formed by collapse, while suffosion dolines develop by suffosion [10, 12].

The development of dropout dolines is also contributed by earthquakes [13, 14]. Earthquake waves may trigger rock failure and thus, the collapse of rocks that became looser.

Irreversible and reversible (fluctuating) water level changes occur to the effect of earthquakes [15], and they also affect collapses.

#### **2. Description of the area**

The Sunja River Valley is situated between two larger geographical/landscape units of Sisak Posavina and Banovina (Banija) (**Figure 1**) and belongs to the Petrinja-Sunja hilly terrain. Permanent surface watercourses constitute a highly developed parallel drainage network with significant amounts of surface water. The Sunja River is a right-bank tributary of the Sava River and is a part of the Sava River Basin (Danube River Basin).

*The Impact of Earthquakes on Dropout Doline (Cover Collapse Sinkhole) Development: A Case… DOI: http://dx.doi.org/10.5772/intechopen.108277*

The Sunja River Valley is a covered karst where subsidence dolines occur, but they are also widespread between Hrvatska Kostajnica and Petrinja. Following the earthquakes of 2020–2021, further subsidence (dropout) dolines developed within a relatively short time in relatively great numbers [16, 17]. The altitude of the bearing surface is 175–190 m. The established level of the groundwater is 177–183 m and

#### **Figure 2.**

*Hydrological map of the environs of the dolines ([18], modified). 1. Gravel, sandy clay, clay, coarse-grained stream load and boulder, 2. Badenian limestone, 3. Sandstone, conglomerate, marl, calcareous and clayey marl, sand, clay, 4. Detected geological boundary, 5. Assumed geological boundary, 6. Erosional-discordance boundary, 7. Detected fracture, 8. Assumed fracture, 9. Fracture series, 10. Doline, 11. Direction of groundwater flow, 12. Water level according to the state of 11 November 2005, 13. Fluvial measurement site, 14. Geological cross section (A-B) and the cross section of its groundwater level (C-D), and geophysical cross section (E-F) 15. Occupied permanent spring with permanent water, 16. Occupied spring with permanent water, for local use, 17. Spring of permanent water with low discharge, 18. Intermittent spring with higher discharge, 19. Intermittent spring, 20. Drilling, 21. Planned drilling, 22. Drilled well, 23. Dug well, 24. Subsurface water level.*

generally flows in S-E and N-E direction towards the springs Pašino vrelo, Bojanića vrelo and Davidovića vrelo (**Figure 2**). These springs create lakes in inactive dolines. The spring water of Pašino vrelo originates from Quaternary sediments in 70%, from Badenian limestone in 30%, while at Bojanića vrelo this proportion is 66% and 34%, at Davidovića vrelo 90% and 10% and at PBV-3 well 56% and 44% [18, 19].

The appearance of the springs is the consequence of a pronounced diagonal fault in the Sunja River Valley, which brought into direct contact permeable Upper Badenian predominantly carbonate rocks especially lithothamnium limestones (M4) and poorly permeable to impermeable Pannonian marls (M6) [18, 19]. Also, the area consists of Plio-Pleistocene age (Pl-Q) constituted by clay, clayey sand and pebbles and the alluvium (al- Q2) constituted by pebbles, clay, clayey pebbles, older rocks, deluvial and proluvial sediments (**Figures 2** and **3**) [18, 20].

Geological structures in the investigated area extend mostly in NNW-SSE or NW-SE directions and follow the so-called 'Dinaric'strike (NW-SE), with predominantly dipslip movements which are tectonically disturbed by the intersection of longitudinal NW-SE right-lateral and transverse NE-SW left-lateral faults of different size and importance along the transitional contact zone of the Dinarides and the Pannonian Basin [8]. The most important fault zone in the area is an active fault zone Pokupsko-Banja Luka in the Dinaridic ophiolite zone, Sava zone. According to [21], the strongest earthquake in the Kupa Valley M = 5.8 was recorded in the year 1909 (Mohorovičić discontinuity). The first strong earthquake in the wider Petrinja area was recorded on 28 December 2020. The day after, the second and third ones had magnitudes of M = 5 and M = 6.2 (h = 10 km). These earthquakes occurred on the Hrastovički fault, a segment of the Pokupsko Fault Zone stretching from Jastrebarsko through Pokupsko towards Banja Luka. These earthquakes caused the opening of approximately 100 new sinkholes (dolines) in the wider Mečenčani and Borojevići area.

Water inflow features are dolines, which are covered karst features and occur on the floor of the Sunja Valley. Among the dolines, there are distinguished [20] old dolines (they developed preceding the earthquake), more recent dolines (which

#### **Figure 3.**

*Geological cross section (A-B) [18]. 1. Alluvium of the Sunja River: gravel, clay, sand, boulder, 2. Badenian limestone, 3. Eocene, Ottnangian, Sarmatian, Pannonian clastic sediments, 4. Boundary of beds, 5. Erosional discordance boundary, 6. Fracture, 7. Structural-piezometric drilling.*

*The Impact of Earthquakes on Dropout Doline (Cover Collapse Sinkhole) Development: A Case… DOI: http://dx.doi.org/10.5772/intechopen.108277*

developed during the earthquakes) and even newer dolines (which developed after the earthquakes) and primitive dolines and buried dolines (**Figures 4** and **5**). Old features (42 dolines) have gentle slopes, are covered with vegetation and without

**Figure 4.**

*Distribution of dolines (modified from [20]). 1. Buried doline, 2. Primitive doline, 3. New doline, 4. Old doline, 5. Development date (year that is missing from the map: in 2020, the month is 12 and in 2021, the months are January and February), 6. Boundary of the area marked for the calculation of the development rate.*

#### **Figure 5.**

*New (dropout) subsidence dolines: A. non-narrowing dropout doline with lake (photo taken by Željko Grgić in February 2021), B. dropout doline with lake that is narrowing towards its floor from close-up (photo taken in February 2021 by Ronald Goršić): 1. Traces of liquefaction, 2. Limestone block, 3. Thrown out material, 4. Ragged lawn with traces of primary collapse, 5. Secondary collapse, C. narrowing dropout doline with lake from far view (a) and dropout doline without water (b) (photo taken in January 2021 by Marijan Car, Mario Bačić, Josip Terzić). D. Twin-like doline: Partial depression at the right side of the photo which developed either by the newer collapse of the cover cavity or during the secondary collapse of the side slope (photo taken by Sonja Zlatovič) 1, 2. Partial depressions, 3. Secondary collapse. All dolines were survey measured and monitored by Jeronim Moharić from January to September 2021 for the purpose of this study.*

water. During the earthquakes, 82 new dolines developed in an area of 4 km<sup>2</sup> (but doline development also continued after the earthquake activity, and by the beginning of May, their number was 91 and by the beginning of December 2021, they numbered more than 100), which are collapse features with steep slopes (**Figure 5**). Several boreholes reached greater depths in order to get information on the composition of the cover. In case of a dropout doline, the sediments of the cover are organic soil, sandy lean clay and lean clay in 8-metre thickness [22].

Water outflow sites are springs and spring lakes and the depressions containing water (paleodolines) such as the spring Pašino vrelo.

#### **3. Methods**


*The Impact of Earthquakes on Dropout Doline (Cover Collapse Sinkhole) Development: A Case… DOI: http://dx.doi.org/10.5772/intechopen.108277*


The calculation of the ceiling thickness that can collapse without earthquakes for a given cavity diameter was done and led to the formation of a dropout doline at the surface. Using this relatively simple theoretical model, it is possible to determine the critical ceiling thickness at which collapse will necessarily occur.

The model under consideration is shown in **Figure 6**. It depicts a cylindrical cavity of diameter *D* at depth *H* below the surface. The initial assumption of model is that the ceiling of the cavity is not subject to any supporting forces from below, so the

#### **Figure 6.**

*Sketch of a theoretical model of dropout dolines. Explanation: 1. Uncollapsed rock mass, 2. Collapsed rock mass, 3. Cavity, 4. Interface between collapsed and intact rock masses, 5. Vertical stress in the rock, 6. Horizontal stress in the rock, 7. Compressive stress along the discontinuity surface, 8. Shear stress along the discontinuity surface, 9. The gravity of the destroyed rock aggregate, 10. The diameter of the cavity, 11. The thickness of the ceiling, 12. The slope of the discontinuity surface, 13. The infinitely thin layer of rock (imaginary) cut from the destroyed aggregate, 14. The depth below the surface of the selected layer, 15. The radius of the cut rock layer (disc), 16. The thickness of the cut rock layer.*

equilibrium of the rock mass above the cavity is ensured by the cohesion force and internal friction between the particles of the overlying aggregate. If these forces cannot counterbalance the weight of the rock mass, collapse will occur.

Experience from mining and soil mechanics shows that in low-strength rocks such as the alluvium of the Mečenčani and Borojevići areas, collapse occurs along the upward-expanding fracture surfaces that start at the edge of the cavity and bound a conoidal frustum. The rock beds within the conoidal frustum are subjected vertically downwards to the gravitational forces of the Earth's gravity and upwards to the cohesive and frictional forces of the rock particles outside the fracture surface along the conoidal frustum mantle. These forces are defined below. The rock beds within the stump cone are vertically downwardly influenced by the gravitational forces due to the Earth's gravity and upwardly by the cohesive and frictional forces of rock particles outside the fracture surface along the conoidal frustum mantle. The formulae defining these forces are derived below.

The weight of the falling rock (*G*)

$$\mathbf{G} = \rho \mathbf{g} \mathbf{V},\tag{1}$$

where *ρ* is the density of the rock [kg/m<sup>3</sup> ], *g* is the acceleration due to gravity (*g* = 9,81 m/s<sup>2</sup> ), *V* is the volume of the rock mass [m<sup>3</sup> ], which in the case of a conoidal frustum is

$$W = \frac{\pi}{12}H \cdot \left(3D^2 + 3H \cdot D \operatorname{ctg} \, a + H^2 \operatorname{ctg}^2 a\right). \tag{2}$$

*α* is the angle of the fracture surface with respect to the horizontal, which, according to the relationship known from soil mechanics [23]

$$a = 4\mathbb{S}^{\circ} + \frac{\phi}{2},\tag{3}$$

where *ϕ* is the angle of internal friction of the rock.

There are two types of forces acting along the conoidal frustum sheath (i.e. the discontinuity surface). One is proportional to the compressive stress (*σ*) perpendicular to the cone surface, the other proportional to the shear stress (*τ*) parallel to the cone's surface. To determine *σ* and *τ*, let us imagine a prismatic body in the rock, infinitesimally small, bounded by a horizontal surface, a vertical surface and a surface parallel to the discontinuity surface. This body is loaded from above (in the vertical direction) by the pressure (*σv*):

$$
\sigma\_v = \rho \text{gz},
\tag{4}
$$

and sideward (on the vertical side of the imaginary prism), a horizontal stress (*σh*) is applied:

$$
\sigma\_h = \frac{\nu}{1 - \nu} \rho \text{gz},
\tag{5}
$$

where *z* is the depth below-ground surface [m], and *ν* is the Poisson's ratio of the rock. (Poisson's ratio is the relationship between the transverse and longitudinal deformation of rock under uniaxial pressure. Its value varies between 0 and 0.5.) The equilibrium of forces acting on the prism requires that stresses (*σ, τ*) also occur on the *The Impact of Earthquakes on Dropout Doline (Cover Collapse Sinkhole) Development: A Case… DOI: http://dx.doi.org/10.5772/intechopen.108277*

inclined side. These give the stresses along the discontinuity surface. The compressive stress, which is perpendicular to the surface (according to the elementary relationship known from the strength theory), is

$$
\sigma = \sin^2 a \cdot \sigma\_h + \cos^2 a \cdot \sigma\_v. \tag{6}
$$

The shear stress (*τ*) acting parallel to the discontinuity surface at the moment of initiation of collapse consists of the internal friction between the rock particles sliding side by side, which is directly proportional to the sum of the compressive stress on either side of the fracture surface and the specific cohesive force "pulling" the particles together:

$$
\mathfrak{r} = \mathfrak{tg}\,\phi \cdot \sigma + \mathfrak{c}, \tag{7}
$$

where *c* is cohesion [N/m<sup>2</sup> ].

The stresses *σ* and *τ* are specific, that is they give the force per unit area of the envelope of cone. Let *fv* denote the vertical component of the resultant of these forces, which, according to elementary calculations, is

$$f\_v = \cos a \cdot \sigma + \sin a \cdot \text{r.}\tag{8}$$

According to terms (3)-(7) (after appropriate aggregations)

$$f\_{\nu} = \left(\nu \frac{\cos a \sin^2 a + \text{tg}\,\phi \sin^3 a}{1 - \nu} + \cos^3 a + \text{tg}\,\phi \sin a \cos^2 a\right) \rho \text{g}z + \sin a \cdot c. \tag{9}$$

Imagine cutting out a 'ring' of height *dz* from the cone's mantle at depth *z*. Its radius (*r*), is obviously

$$r = \frac{D}{2} + (H - z)\coth a.\tag{10}$$

Mark the area of the cloak of this ring *dA*:

$$dA = \frac{2r\pi}{\sin a} dz.\tag{11}$$

Since *fv* is constant along this ring, the resultant of the vertical component of the compressive and frictional forces acting on the rock bed is

$$dF = f\_v \cdot dA,\tag{12}$$

If we add up the forces acting on all the elementary rings covering the envelope of cone of the rock-face before the collapse, we obtain the resultant *F* of the forces acting upwards (opposite to gravity) on the body, which, taking into account (9), (10) and (11)

$$F = \int\_0^h \pi [D + 2(H - z)\circledast a] \cdot \left[ \left( \nu \frac{\cos a \sin a + \text{tg } \phi \sin^2 a}{1 - \nu} + \frac{\cos^3 a}{\sin a} + \text{tg } \phi \cos^2 a \right) \rho gz + c \right] dx \tag{13}$$

For ease of writing, let us introduce the symbol

$$\kappa = \frac{1}{6} \left( \nu \frac{\cos a \sin a + \text{tg}\,\phi \sin^2 a}{1 - \nu} + \frac{\cos^3 a}{\sin a} + \text{tg}\,\phi \cos^2 a \right). \tag{14}$$

Performing the integration formulated in (13)

$$F = \pi\kappa (2H^3 \operatorname{ctg} a + 3DH^2)\rho \mathbf{g} + \pi (H^2 \operatorname{ctg} a + DH)c. \tag{15}$$

Eq. (15) gives the force that ensures the equilibrium of the cavity ceiling. It can be seen that this has a linear (first-order) relationship with cohesion, that is the more cohesive, more frictional rocks have a greater 'holding power' at the same depth than rocks with less cohesion or less internal friction. There is also a linear relationship between the retention force and the density of the rock.

The formation of a doline starts when the balance of forces is upset for some reason, that is when

$$F \le G,\tag{16}$$

The critical (*H*) cavity depth at which collapse can occur is given by *F=G*. Expressions (1) and (14) yield a second-degree equation for *H*:

$$\left(2\text{4\kappa ctg }a - \text{ctg}^2 a\right)\rho \text{g}H^2 + \left[12\sigma\_c \text{tg}\,\phi \text{ ctg}\,a + (3\text{6\kappa}D - 3D \text{ ctg}\,a)\rho \text{g}\right]H + \left(12\text{Dc} - 3D^2 \rho \text{g}\right) = 0\tag{17}$$

A numerical analysis of (17) (with realistic data for the Mečenčani and Borojevići area) shows that the first term is several orders of magnitude smaller than the second and third terms and is therefore negligible. The numerical analysis also shows that the first member of the second bracket term is negligible next to the second and third. The solution to the simplified equation for *H* (taking into account the expression (14) for *κ*) is given by:

$$H = \frac{1 - \nu}{2\nu} \frac{D - \frac{4\epsilon}{\rho \mathbf{g}}}{\sin a (\cos a + \mathbf{t} \mathbf{g} \,\phi \sin a)}. \tag{18}$$

*H* is the critical ceiling thickness (expression 18) at which the roof covering is 'just' not yet collapsing. However, at this ceiling thickness, the cavity is already unstable, that is the slightest disturbance (earthquake, rock scour, soil moisture variation, etc.) will cause immediate collapse.

#### **4. Results**

The pattern of doline distribution shows their tectonic determination. The fault zone detected in the bedrock [3, 17, 19] directed the karst water flow, the transportation of dissolved material and thus, the direction of cavity formation there. However, the cavities that formed in this way determined the site of the material loss of the cover and the zone in which it took place. The geoelectrical resistance values of the profiles indicate the extent of fracturing and cavity formation of the bedrock below the zone of the dolines.

*The Impact of Earthquakes on Dropout Doline (Cover Collapse Sinkhole) Development: A Case… DOI: http://dx.doi.org/10.5772/intechopen.108277*

Based on the earthquake and precipitation data of the area, the following can be established:


Taking into consideration the observed dates of doline development (**Figure 4**, available only for some dolines), the earthquakes may have triggered immediately (on the day of the earthquake) and delayed doline development (days after the earthquakes). An example for immediate doline formation was what happened on 29 December 2020 when two dolines were formed and 30 December 2020 and 31 December 2020 when one doline developed each day. A similar situation was observable on the days of earthquakes of 2 January 2021, 4 January 2021, 5 January 2021 and on 6 January 2021 when one doline was formed each day. On the contrary, a delayed doline formation took place on 23 January 2021 when two dolines developed; however, no earthquake occurred on that day, but they happened only on the days of 7 January 2021, 9 January 2021, 10 January 2021 and 15 January 2021. No doline development coincides with the day of significant amount of rainfall. However, on 29 December 2020 when two dolines developed, the quantity of precipitation was significant. The relation between rainfall (infiltration) and doline development is less obvious. On 29 December 2020, two dolines developed at rainfall of 6 mm, while one doline was formed at precipitation of 12 mm on 30 December 2020. It is more probable that precipitation has a delayed effect. On 26 December 2020 when precipitation was 31 mm, it must have ensured favourable conditions for the development of the above-mentioned dolines on the 29th and 30th of December. The role of precipitation in doline development is indicated by the fact that dolines continued to take shape at the end of 2021. Thus, in December 2021, four new dolines developed probably to the effect of autumn rainfalls and also after a stronger M 3.2 earthquake.

The distribution of depressions is of banded (at some places, of linear) pattern. The band has a NW-SE direction and a length of 3100 m (outside that band only six dolines occur which developed during the earthquakes). Perpendicular to this direction, the width of the band is maximum 800–850 m (**Figure 4**). The band of dolines intersects the Sunja River in NW, and only some new dolines occur on its left bank. The width of the band covered by old dolines is larger; thus, the appearance of new dolines is more concentrated than that of the old dolines. The occurrence of new dolines at the north-western and south-eastern ends of the band (mainly at the latter) is slightly dispersed. However, older dolines frequently group along this arcuate band. The occurrence of dolines is combined in the band: there are sections made up of only new and only old depressions, but mixed sections are also present. Depressions outside this arcuate band are in groups of an irregular pattern. These may be constituted by homogeneous depressions (made up of only new or only old dolines), but by depressions that developed at various time. The direction of the band and the tectonic structure as well as the concordance of the directions of some fractures (faults) refers to the fact that doline development is controlled by tectonics. The area bearing the dolines (this is the floor of the Sunja River Valley) is of a very low inclination. It is 1.29° along the profile perpendicular to the Sunja River (**Figure 2**).

The development rate of new dolines is extremely high. If we consider the area with dolines and presume a period from 28 December to the end of March (doline development is probably shorter than this: in March no doline development was likely to the direct effect of the earthquakes, but the last earthquake, as we observed, took place on 3 March) calculating with 93 days, in the bearing area of 0.0113 km<sup>2</sup> , the average rate of doline development is 27.33 dolines per month. Calculating this to 1 km<sup>2</sup> and 1 month, the rate is 2418.70 doline/month or calculating with the actual area to 1 year, it is 327.96 doline/1 km<sup>2</sup> , calculating to 1 km<sup>2</sup> and 1 year, 201.56 doline develops, and to 1000 m<sup>2</sup> and 1 month the development rate is 29.02 doline/month. Calculating with the actual area, the development rate is 0.91 doline per day. Thus, comparing with the development rates of literary data, in the area of Mečenčani and Borojevići, the same number or more dolines per day were formed during a day than in other karst areas during a year. This high rate of doline development can only be explained by the impact of the series of earthquakes.

The morphological and hydrological characteristics of the dolines are the following:


The depth of new depressions does not exceed 5 m predominantly; thus, they were formed in the unconsolidated superficial deposit (al-Q2). Therefore, they are subsidence dolines. Only one doline has a depth of 12 m. However, it refers to the fact that the collapse triggering its development also spreads over the Badenian limestone. The cavity of the Badenian limestone must have been without water (or partly without water) at least during the time of the earthquake which triggered its collapse since the collapse of the cavity is only possible if the ceiling of the cavity is not supported. However, the development of some depressions (those with a depth close to 5 m) may *The Impact of Earthquakes on Dropout Doline (Cover Collapse Sinkhole) Development: A Case… DOI: http://dx.doi.org/10.5772/intechopen.108277*

have started at the Badenian limestone since the depth of the doline was reduced by the collapsed material.

Among the new ones, there are depressions that have been dry since their development, there may be constantly wet depressions (there is a lake with permanent water on their floor) and there are depressions that are intermittently wet (a lake appears on their floor intermittently). The altitude difference between the water level of the doline with permanent water and the groundwater table is 4.5 m, and the distance between the distance between the doline with lake and the well is 80 m. The similarity in water levels indicates that the water of dolines with lake partly originates from groundwater or the dolines developed from cavities that had at least partly been below the groundwater table. (The reason for the difference of the above water levels may be that the collapse of the cavity increased the water level in the doline and the water level of the lake follows the subsidence of the groundwater level with delay.) The water-level fluctuates in both the permanent and intermittent lakes.

If the altitude of the floor of new dolines is almost the same as the altitude of the cavities from which they developed, then the extent to which the dolines are filled with water or its lack indicates the position of former cavities as compared to the groundwater level. As compared to the groundwater level, the dolines and their cavities may be and may have been of the following types:


Dolines of water outflow sites have steep slopes and permanent lakes. These features are not active (paleodolines) such as the spring Pašino vrelo. They are situated at the emergence sites of groundwater and karstwater. Based on the composition of their water, they are drainage sites for them. According to [19], groundwater from limestones and water from Quaternary deposits are mixed. The groundwater is the water of the al-Q2 beds; in the NE, it follows the Sunja River and is situated at a depth of 1–2 m. Its position as compared to the surface SW from Mečenčani is unknown. However, it can be observed that the altitude difference between the surface and the water level increases farther from the Sunja River (**Figure 7**). The different proportion of the two waters at various drainage sites refers to the existence of two kinds of flow

**Figure 7.** *Position of groundwater. Legend: 1. surface, 2. groundwater level.*

systems, but also to their partial or widespread relation. Their discharge fluctuation was between 24.4 and 12.4 l/s at the spring of Pašino vrelo. The fluctuation of the groundwater discharge and also the fluctuation of its water level [18, 19] refer to the presence of high and low water levels.

The size of dolines that developed at water inflow sites at various times is different. Older dolines are of a greater diameter and a smaller depth (the average diameter is 5.1 m, the average depth is 0.47 m), newer dolines have a smaller diameter, but they are deeper (their average diameter is 3.7 m, and their average depth is 1.5 m). These values are more striking if we only consider the sizes of those with a shape index below 0.3 in case of old dolines (**Figure 8**): here, the average diameter is 5.73 m, and the average depth is 0.34 m. The latter are the smallest, with a diameter of only 1–2 m, and having steep slopes, they are passage-like and aligned with a collapse yard of small depth. The shape of 34 old dolines is in the 0.0–0.1, 0.1–0.2, 0.2–0.3 classes and only eight dolines have a shape larger than class 0.3 (**Figure 8**). The shape of 63 new dolines is in classes larger than 0.3. Thus, at old dolines, a larger diameter belongs to a given depth, while in case of young dolines, a smaller diameter belongs to a given

**Figure 8.** *Shape distribution of old and new dolines and the functions of shape distributions.*

*The Impact of Earthquakes on Dropout Doline (Cover Collapse Sinkhole) Development: A Case… DOI: http://dx.doi.org/10.5772/intechopen.108277*

depth. The class distribution of the two doline groups is also different. The distribution of old dolines can be described with a decreasing polynomial function, while that of new dolines can be described with a polynomial function with two maximum values (**Figure 8**).

New dolines can be referred into two groups. At one of them, the depth is dominant (the depth is larger than the width), at the other, the width is dominant (the width is larger than the depth) since there is a functional relationship between the width and depth at the members (at every depression) of the group (**Figure 9**). If the depressions developed from cover cavities, the former cavities can also be put into two types: at one of them, the vertical expansion is dominant, while at the other, the horizontal expansion is dominant. Therefore, the cavities that were situated at a similar depth as compared to the surface were mainly either vertical or horizontal. During the collapse, the cavities preserved their horizontal and vertical dimensions and they were only modified to a small degree. The small number of vertical depressions (there are four) and the large number of depressions with large width prove that great width favours the denudation of the cavity ceiling.

In case of depressions whose depth exceeds their diameter, there is a strong relation between the width and the depth at the width and depth function (**Figure 9**). Depressions becoming wider and wider are more and more expanded vertically. Thus, former wider cavities were probably more and more expanded downwards. At depressions where their width exceeds their depth, in case of those with a small width, depths are less diverse. At depressions whose width is larger, the diversity of their depth increases. Therefore, in the latter case, the cavities from which the depressions

#### **Figure 9.**

*Functional relationship between the width and depth of new (dropout) dolines. 1. Dolines whose depth exceeds their width, 2. Dolines whose width is larger than their depth, 3. Function fitted to deeper dolines, 4. Function fitted to wider dolines.*


**Table 1.**

*Critical (on a limit of stability) ceiling thickness (*H*) as a function of cavity size (*D*) for different cohesion (*c*) and internal friction (*ϕ*).* *The Impact of Earthquakes on Dropout Doline (Cover Collapse Sinkhole) Development: A Case… DOI: http://dx.doi.org/10.5772/intechopen.108277*

developed had an extremely diverse vertical cavity size during their widening. This is possible since the floors of some cavities from which the depressions developed had a deeper and deeper position in case of similar thicknesses of the cavity ceilings.

According to soil mechanics studies in the Mečenčani and Borojevići area (Tomac et al. 2021c), the alluvium has an angle of internal friction of *ϕ* ≈ 24–28°, a specific cohesive strength *c* ≈ 20–50 kPa, a density of *ρ* ≈ 1800–2000 kg/m<sup>3</sup> , a Poisson's ratio *ν* ≈ 0.25–0.33 and a gravity acceleration *g* = 9.81 m/s<sup>2</sup> (Tomac et al. 2021c). Based on the mean value of these (*ϕ* = 25°, *c* = 20 kPa, *ρ* = 1800 kg/m<sup>3</sup> , *ν* = 0.33), the ceiling may stop above a cavity of diameter *D* = 10 m at a thickness of *H* = 6.7 m.

To illustrate the relationship (18), the evolution of the ceiling thicknesses at the stability limit as a function of the size of the cavity before collapse for different realistic values of cohesion *(c)* and internal friction angle (*ϕ*) is presented in **Table 1**. (*ρ* = 2000 kg/m<sup>3</sup> , *ν* = 0.33) It can be seen that *H* is very sensitive to both *ϕ* and *c.* During earthquakes, these factors can undergo significant changes, which can result in previously quiescent cavity ceilings losing their stability and collapsing.

#### **5. Discussion**

The morphology and shape of the dolines indicate that the dolines of the area between Mečenčani and Borojevići originated by suffosion, collapse and probably compaction.

Old dolines with a shape smaller than 0.3 have a small depth because of limited material transport (suffosion and/or compaction), a large diameter (widespread material transport can take place) and gentle slopes (as a result of slow subsidence and slope denudation). New dolines with a shape larger than 0.3 are deeper because of collapses. Therefore, old dolines with a shape index smaller than 0.3 are suffosion dolines, while new dolines with a shape larger than this are dropout dolines.

Among old dolines occur dolines with a shape larger than 0.3 (**Figure 8**). Therefore, old dolines may also have developed by collapse, but their side slopes became gentle subsequently (pseudo-suffosion doline). Among new dolines, there occur dolines with a shape smaller than 0.3. This may refer to the fact that in case of new dolines of this type, collapse already stopped at an early phase or was followed by subsidence.

In the case of some new (dropout) dolines, the reason for shape increase is not width decrease, but ever greater depth. Presuming a cavity ceiling of small thickness (which is a precondition of collapse), this is possible if the former cavities were increasingly expanded downwards. Thus, cavities develop in the zone of groundwater fluctuation; in addition, in case of those with a shape larger than 1, there was a cavity which was originally in the bedrock; then, its development spreads onto the cover.

The relatively large number of old (suffosion) doline refers to the fact that in this area the tendency to doline development might have been significant earlier as well. This can probably be explained by the low inclination of the bearing terrain in addition to cavity formation that is discussed below. In our case, surface inclination on the valley floor, as already mentioned, is low. This can partly be traced back to the fact that on such terrains, the infiltration proportion of precipitation is large; on the other hand, as a result of its low inclination, the groundwater level is in wide expansion close to the surface (and the karstwater level is also close to the bedrock surface). The water level being close to the surface favours widespread subsurface cavity formation, which enables widespread material transport from the cover. However, collapse

tendency is high because the cover cavities are close to the surface. The cavity formation in the bedrock favours collapses and the reception of the material transported from the cover.

Regarding the water outflow sites of the area, it can be established that the water level oscillates at Pašino vrelo whose degree alternates between 0.46 and 0.52 m since the discharge of its spring fluctuates between 22.8 and 12.8 l/min [18, 19]. This indicates the fluctuation of groundwater level, that is high groundwater level and low groundwater level. Their range between dry and wet seasons is 2 m [16].

Since the groundwater of the cover and the karstwater of the Badenian limestone are present to various proportions at springs Pašino vrelo, at Davidovića vrelo [18, 19], the two kinds of water are separated from each other at least temporarily, and they are primarily mixed at the water outflow sites. The water of the Badenian limestone functions as artesian water in wet seasons [16]. Therefore, the groundwater and the karstwater flow circulate individually.

A two-storey cavity system developed as a result of the two water bodies and their circulation and the fluctuation of groundwater level. The lower level is in the Badenian limestone, while the upper is in the cover.

The shape and size of cover cavities were different. This is proved by the fact that the depth and shape of dolines differ for the dolines (it is presumed that during their development, the dolines preserve the shape and size of former cavities). The shape of the former cavity can be determined from the shape of the depression. The shape, but also the size and position of the cavity point to a collapse tendency. The wider the cavity (the smaller shape it has) and the closer to the surface, the greater its collapse tendency.

Cover cavities are of various positions as compared to the changing groundwater level, since among the developed dropout dolines there occur dry dolines, dolines with permanent lakes and dolines with intermittent lakes. (According to the observation of 05 June 2021, some dolines lost their water.) Among former cavities there were some which were situated in their whole expansion above the high groundwater level and there were others which were above the low groundwater level in their whole expansion. Cavities, whose lower part was below the groundwater table developed by the partial collapse of the bedrock and the cover above it since the permanent presence of groundwater, hindered suffosion.

The altitude of the groundwater level affects the size and development of the cavities and the way of material transport out of the cavities. The difference between the altitude of the groundwater level and the altitude of the surface increases in SW direction from the Sunja River (**Figure 7**). However, the extent of the fluctuation of high and low groundwater levels also increases in the same direction since the water supply from the Sunja River into the groundwater and thus, its affect raising the water level, diminishes moving farther from the river. In case of low water level, its water level reducing effect is limited since the water level cannot reach below the water level of the river.

Consequently, the cavities of the cover may have developed in the following manner:

• By suffosion, cavities in the environs of which the cover is at least temporarily without groundwater and it is coarse-grained. Suffosion is possible because a passage develops above the cavity of the bedrock that reaches above the karstwater level and this transports the material of the cover into the karstic cavity.

*The Impact of Earthquakes on Dropout Doline (Cover Collapse Sinkhole) Development: A Case… DOI: http://dx.doi.org/10.5772/intechopen.108277*


There are a large number of cavities in the upper level of the Badenian limestone according to measurements [18]. The cavities are situated both above and below the water level. Evidence for the existence of cavities above the water level is the already mentioned 12-m-deep-doline. Since it could have developed in a way that it continues in the bedrock, thus, the bedrock cavity had to collapse as well. It is only possible if the cavity reached above the water level.

The material flow of the two-level cavity system may be the following:


There is a higher probability of old (suffosion) doline development at sites where:


New (dropout) dolines may develop in the following ways.


The collapse of cover cavities may take place only during the failure of the ceiling in case of cavities above the high groundwater level (direct earthquake effect) and when it was contributed by other factors too (indirect earthquake effect). The collapses that triggered the development of dolines in the environs of Mečenčani and Borojevići may have been affected by several factors too if the ceiling did not directly collapse because of the failure such as the partial liquefaction of the ceiling material may have happened at cover cavities situated below the high groundwater level. Traces of liquefaction can be seen in **Figure 5B**. Concomitant phenomena, material shots can be recognized in the environs of some dolines (**Figure 5B** and **C**).

In the area of Mečenčani and Borojevići, after the series of earthquakes, the groundwater level became (and stayed) somewhere 30–50 cm lower than earlier in the wells (conversations with locals in 2021). This was probably caused by the irreversible displacement of karstic blocks and the resulting karstwater subsidence spreads onto the groundwater too and/or because of the decrease of cover compaction the groundwater level sank in an irreversible way. The groundwater subsidence reduced the support of cavity ceilings (doline development is without support), which launched a process that is well known in karst literature [24]: the collapse of the ceiling weakened by liquefaction. The thinning out of the ceilings was also enabled by the fact that parts became separated from there, partly because the infiltration of precipitation increased pore water pressure and partly because of liquefaction. However, the infiltrating precipitation also caused an increase in burden, which resulted in the warping of the ceilings. Local people reported on surface subsidence preceding doline formation [16]. This may also have enhanced the separation of parts from the ceilings. To the effect of repeating shocks, the chance of the collapse of thicker ceilings may also increase. Cavities with thin ceilings probably collapsed to the immediate effect of earthquakes, while the cavities with thicker ceilings collapsed due to the delayed effect from the cumulative impact of the shocks.

Dolines, the ceilings of whose cavities had been above the high groundwater level, may have developed during the failure of the cover due to direct earthquake effect by collapse. A favourable condition for this, as already mentioned, is that the cover thickness is less than 1–2 m, and it is not greater than some tens of centimetres. Similarly, the collapse of bedrock cavities was caused by the direct effect of earthquakes.

The relationship (18) highlights several factors that are important for the formation of dropout dolines and even allows a quantitative analysis of these factors.

It can be seen that there is a linear relationship between the thickness (*H*) of the still stable cavity ceiling and the diameter (*D*) of the cavity: the main thickness of a

*The Impact of Earthquakes on Dropout Doline (Cover Collapse Sinkhole) Development: A Case… DOI: http://dx.doi.org/10.5772/intechopen.108277*

cavity with a stable ceiling (i.e. not yet collapsing) is proportionally larger. It is noteworthy that according to this model, there is a minimum cavity diameter below which the ceiling of smaller cavities (due to their own weight) does not collapse. If

$$D < \frac{4c}{\rho \text{g}},\tag{19}$$

then the numerator of (19) becomes negative, which means that no matter how thin the cavity ceiling is, it is held together by cohesion. (Of course, in this extreme case, other external factors neglected in the present model may still trigger collapse.)

From (18), it can be seen that if the cohesion (*c*) of the rock is larger, then *H* is smaller. This is consistent with the experience that more consolidated rocks do not collapse even with thinner ceilings. It also turns out that if the Poisson's ratio (*ν*) of the rock is smaller, a larger cover thickness is required to avoid collapse. The reason for this is because, with a smaller Poisson's ratio, the primary vertical rock pressure (*σv*) results in a smaller horizontal pressure (*σh*) (see relation (4)), which results in a smaller frictional force balancing the cover assembly.

The relationship (18) gives an account of all the effects of earthquakes that lead to the formation of dropout sinkholes. The most important of these (which also play a central role in building damage) are the accelerations induced by the quake waves. These can be calculated from seismology (using the distance from the epicentre or hypocentre of the quake and the magnitude of the quake) and are added to the gravity acceleration. For example, for a 5-magnitude quake, this can increase the value of *g* by up to 30%. Since in expression (18) *g* is contained in the denominator of a subtraction term, and an increase in *g* reduces the subtraction term and therefore increases the value of *H*. In practice, this means that a cavity whose ceiling thickness still ensures stability under 'normal' conditions will no longer comply and collapse due to the increased critical principal strain during an earthquake.

Similarly, the reduction in internal friction and cohesion due to the shear stress of the earthquake increases the value of *H*: a reduction in both *c* and *ϕ* leads to an increase in the value of *H* according to (18); that is, the effective cover thickness becomes less than the minimum value required to avoid collapse.

#### **6. Conclusions**

The karst features of Mečenčani area were classified, the morphology of dolines was described and dolines were put into groups based on their position as compared to the groundwater level. The coincidence of the development of dolines and the time of earthquakes was described. The shape distribution of old and new dolines was also studied.

Statements related to doline development are as follows:


develop by the inheritance of the bedrock cavity onto the surface through the cover.


The study of the dolines of the Mečenčani area draws attention to the fact that doline development is influenced by several factors. Some factors are only specific of this area (the closeness of the karstwater level to the bedrock), while others are generally valid (groundwater, earthquake). The significance of our study is that the role of various factors in doline development was given in quantities. Thus, a better prognosis for doline development can be given on various covered karsts of the Earth and human made structures can be planned more safely.

#### **Acknowledgements**

Many thanks for overall support to Mr. Gordan Hanžek (State Secretary of the Central State Office for Reconstruction and Housing of the Republic of Croatia) and the Ministry of Physical Planning, Construction and State Property. We are very grateful for the help with collection data on the field to Mr. Jeronim Moharić with colleagues (Geo Gauss d.o.o. from Čakovec), Borna Gašpert and Mirko Stanković. The English language was edited by Tatjana Jauk (linguist at Hrvatske vode) and by Dénes Lóczy (University of Pécs), and we thank them a lot for her kind help.

*The Impact of Earthquakes on Dropout Doline (Cover Collapse Sinkhole) Development: A Case… DOI: http://dx.doi.org/10.5772/intechopen.108277*

#### **Conflict of interest**

There is no conflict of interest.

### **Author details**

Márton Veress<sup>1</sup> \*, Natalija Matić<sup>2</sup> , Zoltán Mitre<sup>3</sup> and Gábor Szunyogh<sup>4</sup>

1 Department of Geography, Eötvös Lóránd University, Szombathely, Hungary


© 2022 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

### **References**

[1] Miranda E, Brzev S, Bijelić N, Arbanas Ž, Bartolac M, Bijelić N, et al. Petrinja Croatia. December 29, 2020 Mw 6.4 Earthquake, Joint reconnaissance report, LFE Program of the Earthquake Engineering Research Institute (EERI), Structural Extreme Events Reconnaissance (StEER) Network, ETH Z¨rich, PRJ-2959, 2021. p. 206. DOI: 10.3929/ethz-b-000465058

[2] Tondi E, Blumetti AM, Čičak M, Di Manna P, Galli P, Invernizzi C, et al. Conjugate coseismic surface faulting related with the 29 December 2020. Scientific Reports. 2021;**11**:9150

[3] Bočić N. Structural-geomorphological aspects of the Petrinja earthquake M6.2 (Croatia) preliminary considerations. Hrvatski geografski glasnik. 2021;**83**(1): 5-24. DOI: 10.21861/HGG.2021.83.01.01

[4] Pollak D, Gulam V, Novosel T, Avanić R, Tomljenović B, Hećej N, et al. The preliminary inventory of coseismic ground failures related to December 2020 – January 2021 Petrinja earthquake series. Geologia Croatica. 2021;**74**(1): 189-208. DOI: 10.4154/gc.2021.08

[5] Karimzadeh S, Matsuoka M. A preliminary damage assessment using dual path synthetic aperture radar analysis for the M 6.4 Petrinja earthquake (2021). Remote Sensing. 2021;**13**(12):2267. DOI: 10.3390/rs13122267

[6] Boncio P, Amoroso S, Atanackov J, Baize S, Barbača J et al., Surface faulting during the 29 December 2020 Mw 6.4 Petrinja earthquake (Croatia), EGU General Assembly, EGU21-16575, vEGU21:Gather online. 2021. DOI: 10.5194/egusphere-egu21-16575

[7] Xiong W, Yu P, Chen W, Liu G, Zhao B, Nie Z, et al. Mw 6.4 Petrinja earthquake: a dextral event with large coseismic slip highlights a complex fault system in northwestern Croatia. Geophysical Journal International 2021. 2020;**228**(3):1935-1945. DOI: 10.1093/ gji/ggab440

[8] Markušić S, Stanko D, Penava D, Ivančić I, Bjelotomić Oršulić O, Korbar T, et al. Destructive M6.2 Petrinja earthquake (Croatia) in 2020— Preliminary multidisciplinary research. Remote Sensing. 2021;**13**:1095. DOI: 10.3390/rs13061095

[9] Tomac I, Zlatović S, Athanasopoulos-Zekkos A, Bleiziffer J, Domitrović D, Frangen T, et al. Geotechnical Reconnaissance and Engineering Effects of the December 29, 2020, M6.4 Petrinja, Croatia Earthquake, and Associated S eismic Sequence, Report GEER-072. USA: GEER Association; 2021a. pp. 5-38. DOI: 10.18118/G63T0S

[10] Williams PW, Dolines. In: Gunn J, editor. Encyclopedia of Caves and Karst Science. Routledge, New York, London: Fitzroy Dearborn; 2004. pp. 304-310

[11] Veress M, Covered Karst. Book. Dordrecht: Springer, Springer Geology; 2016. p. 536. DOI: 10.1007/978-94-017- 7518-2

[12] Waltham T, Bell F, Culshaw M. Sinkholes and Subsidence. Berlin, Heidelberg, New York: Springer-Verlag; 2005. p. 382

[13] Yuan D. Environmental and engineering problems of karst geology in China. In: Beck BF, Wilson WL, editors. Karst Hydrogeology, Engineering and Environmental Applications. 1987. pp. 1-11

[14] Gutierrez F, Parise M, De Waele J, Jourde H. A review on natural and

*The Impact of Earthquakes on Dropout Doline (Cover Collapse Sinkhole) Development: A Case… DOI: http://dx.doi.org/10.5772/intechopen.108277*

human-induced geohazards and impacts in karst. Earth Sciences Review. 2014; **138**:61-88

[15] Mádl Szőnyi J, Czaune B, Simon Sz, Erőss A, Zsemle F, Pulay E, et al. Hidrogeológia (Hydrogeology), Digitális Tankönyv, Hungary: ELTE; 2013

[16] Tomac I, Vlahović I, Parlov J, Matoš B, Matešić D, Kosović I, et al. In: Tomac I, editor. Cover-collapse sinkholes. Geotechnical Reconnaissance and Engineering Effects of the December 29, 2020, M6.4 Petrinja, Croatia Earthquake, and Associated Seismic Sequence, Geotechnical Extreme Events Reconnaissance Report GEER-072. USA: GEER Association; 2021b. pp. 52-99. DOI: 10.18118/G63T0S

[17] Tomljenović B, Cvetković M, Sečanj M, Palenik D, Vukovski M, Belić N, et al. Geological engineering, hydrogeological and geophysical investigations for the purpose of defining an effective concept of organized reconstruction in earthquake-affected areas. Study on seismically induced effects of Petrinja series of earthquakes 2020-2021 – preliminary identification of risks. Zagreb: University of Zagreb Faculty of Mining, Geology and Petroleum Engineering; 2021

[18] Mraz V, Dolić S, Kolarić J, Marković T, Dolić M, Košćal S. Water survey works "Pašino vrelo". Croatian Geological Survey. 2005;**2005**:12-49

[19] Larva O, Marković T, Mraz V. Hydrodynamic and hydrochemical conditions at the groundwater source "Pašino vrelo". Geologia Croatica. 2010; **63**(3):299-312. DOI: 104154/gc.2010.24

[20] UZFMGPE, Records of koseismic surface phenomena in the zone of the main seismogenic fault and maps of safety corridors in relation to collapsed sinkholes in the settlements of Mechenčani and Borojevići. Book 2. Zagreb: University of Zagreb , Republic of Croatia Ministry of physical planning, construction and state property. 2021. pp. 10-25

[21] Herak M, Herak D, Markušić S. Revision of the earthquake catalogue and seismicity of Croatia, Terra Nova. 1996; **8**:86-94

[22] Tomac I, Salković I, Kovačević-Zelić B, Domitrović D, Vučenović H, Hrženjak P. Complementary geotechnical and geophysical investigation works. In: Tomac I, editor. Geotechnical Reconnaissance and Engineering Effects of the December 29, 2020. Report GEER-072, USA: GEER Association; 2021c. pp. 164-196. DOI: 10.18118/G63T0S

[23] Hudson J, Harrison J, Engineering Rock Mechanics: An Introduction to the Principles, Second Impression. London: Elsevier Science Ltd; 2000. p. 114

[24] Jia L, Li L, Meng Y, Wu Y, Pan Z, Yin R, Responses of cover-collapse sinkholes to groundwater changes A case study of early warning of soil cave and sinkhole activity on Datansha Island in Guangzhou, China. Environmental Earth Sciences. 2018;**77**:488. DOI: 10.1007/ s12665-018-7603-3

### *Edited by António Vieira and Resat A. Oyguc*

Applied geomorphology aims to understand the constraints that natural dynamics impose on human activities, as well as societal impacts on geomorphic forms and processes. It is therefore concerned with the analysis and interpretation of landforms resulting from the interaction between anthropic and non-anthropic (so-called natural) processes, using methodologies specific to this scientific area. This book provides a comprehensive overview of applied geomorphology. It includes five chapters that address such topics as geodiversity as a tool for nature conservation, geoheritage and its enhancement in the context of geotourism, piles as structural elements, slope stability, and landslides.

Published in London, UK © 2024 IntechOpen © Marina\_Skoropadskaya / iStock

Current Perspectives on Applied Geomorphology

Current Perspectives on

Applied Geomorphology

*Edited by António Vieira* 

*and Resat A. Oyguc*