**SpringerBriefs in Applied Sciences and Technology Lizhou Wu · Jianting Zhou**

# **Rainfall Infiltration in Unsaturated Soil Slope Failure**

# **SpringerBriefs in Applied Sciences and Technology**

SpringerBriefs present concise summaries of cutting-edge research and practical applications across a wide spectrum of fields. Featuring compact volumes of 50 to 125 pages, the series covers a range of content from professional to academic.

Typical publications can be:


SpringerBriefs are characterized by fast, global electronic dissemination, standard publishing contracts, standardized manuscript preparation and formatting guidelines, and expedited production schedules.

On the one hand, **SpringerBriefs in Applied Sciences and Technology** are devoted to the publication of fundamentals and applications within the different classical engineering disciplines as well as in interdisciplinary fields that recently emerged between these areas. On the other hand, as the boundary separating fundamental research and applied technology is more and more dissolving, this series is particularly open to trans-disciplinary topics between fundamental science and engineering.

Indexed by EI-Compendex, SCOPUS and Springerlink.

Lizhou Wu · Jianting Zhou

# Rainfall Infiltration in Unsaturated Soil Slope Failure

Lizhou Wu Chongqing Jiaotong University Chongqing, China

Jianting Zhou Chongqing Jiaotong University Chongqing, China

ISSN 2191-530X ISSN 2191-5318 (electronic) SpringerBriefs in Applied Sciences and Technology ISBN 978-981-19-9736-5 ISBN 978-981-19-9737-2 (eBook) https://doi.org/10.1007/978-981-19-9737-2

© The Author(s) 2023. This book is an open access publication.

**Open Access** This book is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this book are included in the book's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the book's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

The use of general descriptive names, registered names, trademarks, service marks, etc. in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant protective laws and regulations and therefore free for general use.

The publisher, the authors, and the editors are safe to assume that the advice and information in this book are believed to be true and accurate at the date of publication. Neither the publisher nor the authors or the editors give a warranty, expressed or implied, with respect to the material contained herein or for any errors or omissions that may have been made. The publisher remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This Springer imprint is published by the registered company Springer Nature Singapore Pte Ltd. The registered company address is: 152 Beach Road, #21-01/04 Gateway East, Singapore 189721, Singapore

# **Foreword by Jianbing Peng**

It is a great pleasure to write this Foreword for Dr. Lizhou Wu's book on *Rainfall Infiltration in Unsaturated Soil Slope Failure*. I have ever taken a part in guiding the research and enjoyed our discussions on rainfall-caused landslides.

Rainfall-induced landslides pose hazards to human safety and economic development in the world. The complex geological environments and variable precipitation lead to difficulties in the assessments of slope stabilities. Lizhou's research group has performed a lot of fundamental researches on the scientific problems of rainfall-induced slope failure.

This book by Lizhou Wu and Jianting Zhou bring many creative studies on rainfall seepage and slope failure mechanisms using both numerical techniques and analytical methods. The authors are particularly successful in addressing the analytical and numerical solutions of linear and nonlinear infiltrations related to the investigation of rainfall-induced unsaturated soil slope failures. It proposes creative novel methods for investigating the linear and nonlinear systems of rainfall infiltration inducing unsaturated soil slope failure. The analytical solutions to water infiltration consider vegetation root and coupling effects, which deserve further research. These methods in the book are applicable not only to unsaturated infiltration issues, but also to rainfall-caused landslides.

I am pleased to witness the publication of this book. Students, engineers, and researchers who focus on the unsaturated infiltration and slope stability can benefit a lot from this book. I am confident that this book will serve as an important source of reference for future studies on rainfall-caused landslides.

Jianbing Peng Academician of Chinese Academy of Sciences Chang'an University Xi'an, China

# **Foreword by Huiming Tang**

The book titled *Rainfall Infiltration in Unsaturated Soil Slope Failure* by Lizhou Wu and Jianting Zhou is a welcomed addition that provides an engineering methodology for a challenging problem of interest the world over. Virtually, many countries encounter serious geological hazards including rainfall-triggered landslides.

The book brings primary application areas of engineering which have been individually successful within civil engineering, namely the modeling of rainfall infiltration and rainfall-induced landslides. The authors are well aware of the extensive research that has been undertaken in various parts of the world related to modeling rainfall infiltration. The development of numerical and analytical methods of rainfall seepage and slope stability is particularly well presented in this book. The authors are particularly successful in addressing the many components related to the assessment of rainfall-induced landslides.

When addressing issues related to modeling rainfall-induced landslides, there is a need to understand reliable numerical solutions of the nonlinear seepage equations. The authors have also made a significant contribution in developing numerical methods related to the estimation of unsaturated seepage for investigating rainfallinduced landslides. I believe that this book will also form an important reference book that will be in demand in university and other engineering libraries.

Huiming Tang Chair Professor, China University of Geosciences Wuhan, China

# **Foreword by Limin Zhang**

Rainfall-triggered landslides are a common disaster around the world. Their mechanisms are still not well understood and their evaluating methods not well developed. Rainfall infiltration is particularly difficult to assess when soil slopes are in complex geological environments. Many methods including analytical and numerical solutions have been developed to predict slope instability due to rainfall infiltration. The predicted slope performance may deviate from the reality because of neglecting complex infiltration processes and geoenvironments. Due to the high nonlinearity of the equations governing rainfall infiltration into slopes and the hydraulic property functions, the efficient and reliable solution of the infiltration equations related to rainfall-triggered landslides is a difficult task.

This book by Lizhou Wu and Jianting Zhou combines the complex topics of analytical and numerical modeling of rainfall infiltration and seepage into unsaturated soil slopes along with the analysis of slope stability. Particular attention is given to understanding nonlinear unsaturated infiltration process. The authors have brought the many components related to rainfall-induced landslides in one nutshell. The authors have also made a significant contribution to the development of numerical and analytical methods related to the evaluation of unsaturated infiltration into soil slopes for stability analysis purposes. I am confident that this book will be an important reference to researchers, graduate students, and practicing engineers.

Limin Zhang Chair Professor, Head of Department of Civil and Environmental Engineering The Hong Kong University of Science and Technology Hong Kong, China

# **Preface**

Landslides pose immediate threats to the infrastructures such as buildings, roads, bridges, and geoenvironment in the mountainous areas. Rainfall is the primary trigger of landslides that frequently cause fatalities and large economic loss. Rainfall-induced landslides gain greater attentions as climate becomes more extreme.

A number of complicated mechanisms of rainfall-caused landslides are involved in the analysis of slope stability due to rainfall. Studies on rainfall-induced landslides require good knowledge of not only mechanical properties of the soil, but also hydrologic behavior of the soil governed by the soil seepage properties.

The primary objective of the book clearly presents rainfall infiltration and slope stability analysis methods using analytical and numerical approaches. Analytical solution can consider coupled infiltration and deformation in unsaturated soil slopes considering vegetable root. A series of improved linear or nonlinear iterative methods are developed to solve complex nonlinear infiltration equation, which improves the convergence rate, accuracy, and stability. These methods are applied to simulate rainfall infiltration-induced unsaturated soil slope failures.

The improved numerical methods, nonlinear and linear iterative methods which can be used to address the related unsaturated infiltration problems are also presented. This book is an essential reading for researchers and graduate students who are interested in rainfall infiltration, slope stability, landslides, and geohazards in the fields of civil engineering, engineering geology, and earth science. The book is written to guide professional engineers and practitioners in slope engineering and geohazard management. This book can enhance their understanding of rainfall-induced landslides, help them analyze a specific problem, and prevent landslides and design engineering slopes according to the local soil and climate conditions.

> Lizhou Wu Jianting Zhou Chongqing Jiaotong University Chongqing, China

**Acknowledgements** The authors would like to express the special appreciation and respect to Prof. Jianbing Peng, who gives valuable guidance, constructive suggestions, and continued encouragement throughout this research work. The first author is grateful to his supervisors, Prof. Runqiu Huang and Prof. Limin Zhang, for years of supervision and help. The first author also extends thanks to Prof. Hong Zhang and Prof. Jun Yang for their encouragement and help.

The authors are grateful to Dr. Shuairun Zhu for his contribution to Chaps. 3 and 4, and the editing help provided by Mr. S. H. Li, Mr. Bo He, Mr. Hao Li, and Mr. Tao Ma.

The support of the National Natural Foundation of China (Nos. U20A20314, 41790440, 41790432, and 42277183) and the Natural Science Fund for Distinguished Young Scholars of Chongqing (cstc2020jcyj-jqX0006) are acknowledged.

# **Contents**



# **Chapter 1 Background**

Heavy rainfall in extreme climates often causes natural disasters such as floods, landslides, and debris flows. Rainfall-induced slope instabilities are major geological natural disasters (Glade 1998; Dai et al. 1999; Iverson 2000; Lee and Pradhan 2007; Li et al. 2016a, b; Wu et al. 2020) that can result in considerable loss of life and damage to infrastructure. Extreme events such as storms, which are becoming more severe because of climate change, can trigger fatal landslides. Storm-induced slope failures frequently occur because of rainfall infiltration, particularly in tropical areas (Fourie 1996; Cevasco et al. 2014). Global climate change in many mountainous areas could lead to more severe fluctuations in rainfall, and trigger of soil slope deformations and even slope instability because of the alteration of intensity, frequency, and quantity of rainfall (Dixon and Brook 2007; Jeong et al. 2008). The influence of climate change on rainfall characteristics has the potential to alter the stability of unsaturated soil slopes. Rainfall infiltration causes a decrease in matric suction and an increase in moisture content and hydraulic conductivity in unsaturated soils. The rainfall intensity and duration, initial water table, and hydraulic conductivity are the parameters that significantly affect slope stability (Ng and Shi 1998). An increase in pore-water pressure can reduce the effective stress and thereby weaken the shear strength of slopes. Complex geological environment and human engineering activities are also significant factors of slope instability under rainfall conditions.

Rainfall-induced slope failures have been examined based on experimental modeling, analytical and numerical methods (Ng and Shi 1998; Iverson 2000; Chen et al. 2005; Wu et al. 2009, 2016, 2020; Zhu et al. 2019). Laboratory and field experiments have been carried out to examine the infiltration mechanisms associated with rainfall-induced slope failures (e.g., Lee et al. 2011; Wu et al. 2015, 2017, 2018). Many numerical and analytical studies have investigated the hydraulic responses of slopes to rainfall infiltration and the stability of slopes under such conditions (Iverson 2000; Cai and Ugai 2004; Wiles 2006; Ali et al. 2014; Zhu et al. 2022). Numerical analysis solves for the matric suction or pressure head distribution in a soil slope with

<sup>©</sup> The Author(s) 2023

L. Wu and J. Zhou, *Rainfall Infiltration in Unsaturated Soil Slope Failure*, SpringerBriefs in Applied Sciences and Technology,

varying permeability, and considering different surface conditions of a soil. Many cases demonstrate that change in rainfall patterns may lead to slope failure due to infiltration. The results can provide an indication of the potential influence of climate change on shallow landslides in many mountainous areas (Kim et al. 2012).

Slope stability problems are commonly encountered in engineering projects. Many slope failures are attributed to water infiltration (Cho and Lee 2002; Cai and Ugai 2004). Matric suction is crucial to the stability of soil slopes because dissipation of matric suction leads to decrease in shear strength of unsaturated soils. Slope failure is closely related to the rainfall-induced transient infiltration of slopes in unsaturated soils (Fredlund and Rahardjo 1993; Wu et al. 2020). Shallow landslides are related to periods of intense rainfall. Engineering activities can result in severe geological and environmental issues. Slope failures induced by engineering activities may occur and may progress into landslides. The internal mechanics of slope movements are stress redistribution and the consequent changes in engineering–geological conditions (Marschalko et al. 2012).

Landslide forecasting should take into account rainfall infiltration into soil slopes. The models in predicting the timing and location of landslides are related to dynamic water infiltration in soil slope. The coupled process in an unsaturated soil is of major interest because of its implications for disaster prevention and environmental issues. Analytical approaches have been developed to provide a basic understanding of unsaturated infiltration in terms of the coupling effect (Wu et al. 2020). Meanwhile, numerical approaches provide a powerful tool for solving complex, nonlinear infiltration into unsaturated soils. These numerical models effectively investigate the coupled hydro-mechanical problem involved in unsaturated rooted slope stability issues (Oka et al. 2010; Wu et al. 2020; Zhu et al. 2022).

### **1.1 Rainfall Infiltration Equation**

According to the modified Green–Ampt model (Mein and Larson 1973), the water infiltration process in soils during uniform precipitation can be divided into two stages: a stage controlled by rainfall intensity and a stage controlled by pressure head. The infiltration rate (*f* a) determined by rainfall intensity (*q*r) can be written as:

$$f\_{\mathbf{a}} = q\_{\mathbf{r}} \cos \beta \,\tag{1.1}$$

where *q*r is the rainfall intensity; and β is the slope angle.

The infiltration rate (*f* b) determined by the pressure head can be expressed as:

$$f\_{\rm b} = k\_s \frac{z\_{\rm f} \cos \beta + s\_{\rm f}}{z\_{\rm f}} \tag{1.2}$$

where *k*s represents the permeability coefficient at saturation; *s*f represents the suction head of the wetting front; and *z*f represents the wetting front depth.

If *f* a = *f* b, the water ponding time (*t*p) can be obtained as:

$$t\_{\rm p} = \frac{s\_{\rm f}(\theta\_{\rm s} - \theta\_{\rm i})}{\cos^2 \beta (q\_{\rm r}/k\_{\rm s} - 1)q\_{\rm r}} \tag{1.3}$$

where θs and θi are the saturated moisture content and initial moisture content, respectively.

According to the water balance principle and Darcy law, wetting front movement during rainfall can be calculated as follows:

$$\begin{cases} z\_{\rm f} = [q\_{\rm r} \cos \beta] \frac{t}{(\theta\_{\rm s} - \theta\_{\rm i})}, & t < t\_{\rm p} \\ \frac{dz\_{\rm f}}{dt} = \frac{1}{(\theta\_{\rm s} - \theta\_{\rm i})} \Big[ k\_{\rm s} \frac{z\_{\rm f} \cos \beta + s\_{\rm f}}{z\iota} \Big], & t \ge t\_{\rm p} \end{cases} \tag{1.4}$$

Equation (1.4) is the GA model suitable for slopes (Chen and Young 2006).

A penetration test of saturated sand layers was conducted and found a quantitative relationship between the infiltration velocity of water in the soil and the head loss, namely Darcy law:

$$q = -k\_s \frac{\mathrm{d}H}{\mathrm{d}L} \tag{1.5}$$

where *q* is the flux (discharge per unit area, with units of length per time, m/s), *H* is the total water head, and *L* is the seepage length.

When soil mass is unsaturated, most scholars believe that Darcy law can also be used for analyzing water movement in unsaturated soils. Richards (1931) extended the Darcy law of saturated soil to the unsaturated infiltration and introduced the continuity equation to derive the equation of motion of the unsaturated soil–water flow, namely the Richards' equation. The permeability coefficient (*k*) is expressed as a function of matric suction, and Darcy's law can be expressed as:

$$q = k \nabla H \tag{1.6}$$

where ∇ *H* is a hydraulic gradient, including two components of gravity and suction.

Continuity equation:

$$\frac{\partial \theta}{\partial t} = -\nabla \cdot q \tag{1.7}$$

Then:

$$\nabla \cdot [k \nabla H] = \frac{\partial \theta}{\partial t} \tag{1.8}$$

Substituting Eq. (1.6) into Eq. (1.8), one can obtain:

$$\nabla \cdot \left[ k(h) \nabla (h+z) \right] = \frac{\partial \theta}{\partial t} \tag{1.9}$$

As shown in Fig. 1.1, the Richards' equation governing one-dimensional vertical infiltration in unsaturated soils can be written as:

$$\nabla \cdot \left[ k(h) \nabla h \right] + \frac{\partial k(h)}{\partial z} = \frac{\partial \theta}{\partial t} \tag{1.10}$$

### **1.2 Infiltration Equation for Unsaturated Slopes**

The 2D generalized RE for unsaturated infiltration is expressed as (Ku et al. 2017; Wu et al. 2020):

$$\frac{\partial}{\partial x}\left[K\_x(h)\frac{\partial H}{\partial x}\right] + \frac{\partial}{\partial z}\left[K\_z(h)\frac{\partial H}{\partial z}\right] = C(h)\frac{\partial H}{\partial t} \tag{1.11}$$

where *Kz*(*h*) and *Kx* (*h*) are the permeability coefficients along the vertical direction and lateral direction in unsaturated soils, respectively. To study the groundwater flow of the unsaturated slope (Fig. 1.2), the RE needs to be rotated. Total head can be written as:

$$H = h + E \tag{1.12}$$

and elevation head *E* can be described as:

**Fig. 1.1** 1D rainfall infiltration model

**Fig. 1.2** A slope infiltration model

$$E = x\sin\beta + z\cos\beta\tag{1.13}$$

By substituting Eqs. (1.12) and (1.13) into Eq. (1.11), Eq. (1.11) can be re-expressed as:

$$\frac{\partial}{\partial x}\left[K\_x(h)\left(\frac{\partial h}{\partial x} - \sin \beta\right)\right] + \frac{\partial}{\partial z}\left[K\_z(h)\left(\frac{\partial h}{\partial z} + \cos \beta\right)\right] = C(h)\frac{\partial h}{\partial t} \tag{1.14}$$

According to Iverson's model (2000), the modified RE only considers infiltration in the vertical direction, which can be given by:

$$\frac{\partial}{\partial z} \left[ K\_z(h) \left( \frac{\partial h}{\partial z} + \cos \beta \right) \right] = \frac{\partial \theta}{\partial t} \tag{1.15}$$

### **1.3 Linearized Richards' Equation**

Combined with an exponential model, Eq. (1.15) is linearized. Here, a new parameter *h*<sup>∗</sup> is defined as:

$$h^\* = \mathbf{e}^{ah} - \lambda \tag{1.16}$$

where λ is the key parameter for solving the linearized Richards' equation in the conversion method, which can be defined as a constant λ = e<sup>α</sup>*<sup>h</sup>*d (Tracy 2006; Liu et al. 2015); *h*d is the pressure head value when the soil is dry. The relative permeability coefficient is expressed as:

$$K(h) = \frac{K\_\varepsilon(h)}{k\_s} = \mathbf{e}^{ah} \tag{1.17}$$

Taking the derivative of Eq. (1.16) with respect to *z*, one can obtain:

$$\frac{\partial h^\*}{\partial z} = \alpha \mathbf{e}^{ah} \frac{\partial h}{\partial z} \tag{1.18}$$

Equation (1.18) can be re-stated as follows:

$$\frac{\partial h}{\partial z} = \frac{1}{\alpha} \mathbf{e}^{-\alpha h} \frac{\partial h^\*}{\partial z} \tag{1.19}$$

Furthermore, substituting Eq. (1.16) into Eq. (1.19), one can obtain:

$$K\frac{\partial h}{\partial z} = \mathbf{e}^{ah} \left(\frac{1}{\alpha} \mathbf{e}^{-\alpha h} \frac{\partial h^\*}{\partial z}\right) \tag{1.20}$$

Equation (1.20) can be rewritten as:

$$K \frac{\partial h}{\partial z} = \frac{1}{\alpha} \frac{\partial h^\*}{\partial z} \tag{1.21}$$

Equation (1.17) can be derived from *z* again:

$$\frac{\partial K}{\partial z} = \alpha \mathbf{e}^{\alpha h} \frac{\partial h}{\partial z} \tag{1.22}$$

Substituting Eq. (1.19) into Eq. (1.22), one can have:

$$\frac{\partial K}{\partial z} = \frac{\partial h^\*}{\partial z} \tag{1.23}$$

An exponential model is employed to describe soil moisture (Gardner 1958):

$$\theta(h) = \theta\_{\mathbf{r}} + (\theta\_{\mathbf{s}} - \theta\_{\mathbf{r}}) \mathbf{e}^{\alpha h} \tag{1.24}$$

where θ (*h*) is volumetric water content; θr is the residual volumetric water content.

The derivation of both sides of Eq. (1.24) with respect to time *t* has the following relationship:

$$\frac{\partial \theta}{\partial t} = (\theta\_\mathrm{s} - \theta\_\mathrm{r}) \frac{\partial K}{\partial t} = (\theta\_\mathrm{s} - \theta\_\mathrm{r}) \frac{\partial h^\ast}{\partial t} \tag{1.25}$$

Substituting Eqs. (1.19), (1.23), and (1.24) into Eq. (1.10), the linearized Richards' equation can be obtained:

$$\frac{\partial^2 h^\*}{\partial z^2} + a \frac{\partial h^\*}{\partial z} = c \frac{\partial h^\*}{\partial t} \tag{1.26}$$

where *c* = α(θs − θr)/*k*s.

Equation (1.26) is also expressed as:

$$K\_a \frac{\partial^2 h^\*}{\partial z^2} + K\_\theta \frac{\partial h^\*}{\partial z} = \frac{\partial h^\*}{\partial t} \tag{1.27}$$

where *K*<sup>θ</sup> = *k*s/(θs − θr), *Ka* = *K*<sup>θ</sup> /α.

The finite difference format of the linearized RE (Eq. 1.26) can be expressed as:

$$K\_a \frac{h\_{i+1}^{\*n} - 2h\_i^{\*n} + h\_{i-1}^{\*n}}{\Delta z^2} + K\_\theta \cos \beta \frac{h\_{i+1}^{\*n} - h\_{i-1}^{\*n}}{2\Delta z} = \frac{h\_i^{\*n} - h\_i^{\*n-1}}{\Delta t} \tag{1.28}$$

where *i*, Δ*z*, Δ*t*, and *n* denote the nodal point, grid size, time step, and time level, respectively.

It can be seen from Eqs. (1.26) and (1.27) that the nonlinear partial differential equation (Eq. 1.26) has been transformed into a linear partial differential equation. Once the linear partial differential equation is solved to obtain a numerical solution, the actual pressure head can be written as:

$$h(z,t) = \frac{1}{\alpha} \ln(h^\*(z,t) + \lambda) \tag{1.29}$$

### **1.4 Unsaturated Soil Slope Stability Under Rainfall**

Shear strength is a fundamental material property that is required to address a variety of engineering problems including bearing capacity, slope stability, lateral earth pressure, pavement design, and foundation design. Recently, many researches have focused on the shear strength of unsaturated soils (Fredlund et al. 1996; Lu et al. 2010).

According to Mohr–Coulomb criterion and effective stress, the shear strength of saturated soils can be expressed as:

$$
\pi\_{\mathbb{I}} = c' + \sigma' \tan \varphi' \tag{1.30}
$$

where τf is the shear strength (kPa), *c* is the cohesion (kPa), σn is the normal stress acting on the failure surface (kPa), and ϕ is the angle of internal friction (°). Cohesion and cohesive shear strength are due to chemical bonding between soil particles and surface tension within the water films (Lu and Likos 2006). Frictional shear strength (σn tan ϕ) is owing to internal friction between soil particles that depends on the normal stress acting on the failure surface.

Engineering practices indicate that the shear strength equation of saturated soils can meet the engineering requirements. The shear strength parameters are also influenced by matric suction. With an increase in matric suction, *c* and ϕ increase, which depends on soil texture and structure. Soil shear strength significantly increases with an increase in net normal stress, matric suction, and the parameters of shear strengths.

However, several phases of unsaturated soils make the shear strength equation of saturated soils difficult to apply. Therefore, some studies on the shear strength criteria of unsaturated soils have been carried out. There are main representative shear strength criteria here.

Bishop (1959) developed a shear strength criterion for unsaturated soils:

$$\pi\_{\rm f} = c' + \left[ (\sigma - \mu\_{\rm a})\_{\rm f} + \chi (\mu\_{\rm a} - \mu\_{\rm w})\_{\rm f} \right] \tan \varphi' \tag{1.31}$$

where τf is the shear strength of unsaturated soils; *c*' and ϕ' are the effective cohesion and friction angle, respectively; (*u*a − *u*w) is the matric suction; *u*a is the pore air pressure; *u*w is the pore-water pressure (*h* = *u*w/γ w, γ <sup>w</sup>= ρw*g*); and χ is the function of the degree of saturation.

Based on two stress state variables, the following equation was developed to describe shear strength (Fredlund and Rahardjo 1993):

$$\tau\_{\rm f} = c' + (\sigma - \mu\_{\rm a})\_{\rm f} \tan \varphi' + (\mu\_{\rm a} - \mu\_{\rm w})\_{\rm f} \tan \varphi^b \tag{1.32}$$

where ϕ*<sup>b</sup>*is the internal friction angle due to the distribution of matric suction.

Lu and Likos (2004) proposed a unified form of shear strength equation:

$$\begin{split} \pi\_{\mathfrak{l}} &= c' + \chi\_{\mathfrak{l}} (\sigma - u\_{\mathfrak{a}})\_{\mathfrak{l}} \tan \varphi' + \chi\_{\mathfrak{l}} (u\_{\mathfrak{a}} - u\_{\mathfrak{w}})\_{\mathfrak{l}} \tan \varphi' \\ &= c' + c'' + (\sigma - u\_{\mathfrak{a}})\_{\mathfrak{l}} \tan \varphi' \end{split} \tag{1.33}$$

in which

$$c'' = \chi\_{\rm f}(\mu\_{\rm a} - \mu\_{\rm w})\_{\rm f} \tan \varphi' \tag{1.34}$$

The first two terms in Eq. (1.33), *c*' and *c*'', represent shear strength due to the so-called apparent cohesion in unsaturated soils. In an unsaturated soil, the third term represents frictional shearing resistance provided by the effective normal force at the grain contacts. The apparent cohesion captured by the first two terms includes the classical cohesion *c*' representing shearing resistance arising from interparticle physicochemical forces, and the second term *c*'' describing shearing resistance arising from capillarity effects. The term *c*'' is defined as capillary cohesion hereafter. Physically, capillary cohesion describes the mobilization of suction stress χ(*u*a − *u*w) in terms of shearing resistance. The relationship between capillary cohesion and the maximum suction stress at failure, χf(*u*a − *u*w)f is defined as shear strength also affects the water movement of the soils (Eudoxie et al. 2012).

Slope failure in unsaturated soil regions induced by rainfall is due to shear strength of unsaturated soils (Fredlund and Rahardjo 1993; Lu and Likos 2004; Guzzetti et al. 2008; Muntohar and Liao 2009). Both rainfall characteristics (rainfall intensity and duration) and soil permeability may influence failure mechanism.

The soil slope stability was commonly followed by stability analysis according to the pressure head and/or the stress condition within the soil slope profile. Various techniques were employed to compute factor of safety (*F*s), and the conventional limit equilibrium methods (Alonso et al. 2010). The limit equilibrium approach is mostly effective for slope failure with a small depth compared with their length and breadth. A slope sliding at a depth happens as the driving stress contributing to failure exceeds the anti-slip stress offered by the soil mass strength. Namely, sliding can occur at a particular depth as follows:

$$F\_s(z,t) = \frac{\tan\varphi'}{\tan\beta} + \frac{c'-h(z,t)\chi\_\text{w}\tan\varphi^b}{W\sin\beta\cos\beta} \tag{1.35}$$

where *F*s(*z*, *t*) is the safety factor over depth and time; *W* is the weight of the sliding mass; and γ w represents the unit weights of water.

Equation (1.35) can be re-arranged as:

$$h(z,t) = \frac{c'}{\chi\_{\text{w}} \tan \varphi^b} - \frac{\chi\_{\text{sat}} z \sin \beta \cos \beta}{\chi\_{\text{w}} \tan \varphi^b} \left(F\_{\text{s}} - \frac{\tan \varphi'}{\tan \beta} \right) \tag{1.36}$$

in which, γ sat represents the unit weight of the saturated soil. When *F*s approaches 1, the infinite soil slope reaches a limit state. Based on Eq. (1.36), the limit-state pore-water pressure head can be obtained.

Rainfall-induced landslides may occur in unsaturated soils above the groundwater table, usually with shallow sliding surfaces parallel to the slope surface (Lu and Godt 2008), which involves 2D and 3D problems. However, an infinite slope model is usually used as a simplified model of the 2D or 3D issues with simple geometry and ignores the stress concentration, the practice sometime demonstrates its effectiveness for assessing shallow slope stability (Michalowski 2018).

Slope instabilities are often hydrologically initiated by the advancement of the wetting front alone (Muntohar and Liao 2010), a rise in groundwater level (Asch et al.1999; Montgomery et al. 2009), and positive pore-water pressure on the soil– bedrock boundaries (Baum et al. 2010). The most common mechanism for rainfallinduced landslides occurs when the soil slides on a low-conductivity layer. Rainfall infiltration leads to a rise in the pressure head, resulting in positive pore-water pressures (Iverson 2000; Muntohar and Liao 2010).

Generally, unsaturated soil slope failures happen most frequently during or after rain periods (Wu et al. 2020). The characteristics of water flow, change of porewater pressure, and shear strength of soils are the major parameters related to slope stability analysis involving unsaturated soils that are directly affected by the boundary conditions (i.e., infiltration and evaporation) at the soil–atmosphere interface. The relative importance of soil properties, rainfall intensity, initial water table location, and slope geometry in inducing instability of soil slopes under different rainfall was investigated through a series of studies. Soil properties and rainfall intensity were found to be the primary factors controlling the slope instability due to rainfall, while the initial water table location and slope geometry only played a secondary role (Rahardjo et al. 2007).

The Green–Ampt model is a typical approximate infiltration model. Due to the simplicity and few parameters, the approximate infiltration model has become popular (Grimaldi et al. 2013). The classic GA model is only suitable for infiltration in horizontal soils. Therefore, modified GA models have been developed to describe the water infiltration in layered soils and slopes (e.g., Mein and Larson 1973; Chen and Young 2006; Kale and Sahoo 2011). Some modified infiltration models that account for rainwater redistribution have also been proposed (e.g., Corradini et al. 1997; Dou et al. 2014). These infiltration models have been extended to regional rainfall-runoff models for the hydrological prediction of catchments (Yuan et al. 2019). However, the actual infiltration process is very complicated and affected by many factors such as soil heterogeneity and rainfall conditions, and becomes difficult to be described accurately based on theoretical formulations (Srivastava et al. 2020). These theoretical equations generally tend to overestimate the factor of safety of soil slopes, resulting in slides and geological hazards (Kim et al. 2012). Some intelligent methods have been developed to predict the water infiltration into soils using machine learning techniques (Sihag et al. 2018).

Hydrological responses and slope factor of safety due to rainfall are concerned from a perspective of hydro-mechanical coupling. Coupled and uncoupled hydromechanical behaviors in unsaturated soils have been carried out to characterize the physical responses of unsaturated infiltration (i.e., variation of soil moisture, matric suction, effective stress, shear strength, and slope stability) (Casini 2013). The coupled issues are strongly linked in unsaturated soil slopes due to water infiltration, and the coupled poromechanical model actually examines the behavior and stability of rooted soils subjected to rainfall (Kim et al. 2012). Pressure heads generated in the uncoupled analysis are employed to examine deformation or soil slope stability (Cai and Ugai 2004; Yoo and Jung 2006). The accuracy and computational efficiency of the uncoupled analysis highly depend on the selected time increments (Huang and Lo 2013). The soil hydraulic and mechanical responses are calculated simultaneously in the coupled analysis. The coupled analysis produced a reasonably well defined wetting front and a lower critical *F*s for unsaturated soil slopes. The coupled investigation could produce more accurate assessment of soil slope stability due to water infiltration and demonstrate a better physical representation of water infiltration and stress variation within unsaturated soil slopes.

More and more methods highlight the role of vegetation because of their interception role of the canopy and the root characteristics. Meanwhile, recent studies indicate that vegetation cannot control the rainfall-induced shallow landslide distribution (Emadi-Tafti et al. 2021). Some researches focus on the effect of roots on root–soil composite strength, or saturated hydraulic conductivity (Alessio 2019). The more complex the root architecture is, the stronger the root-composite strength becomes, while the faster the rainfall infiltrates. It has generally been concluded that vegetation roots mechanically and hydrologically affect slope stability. The plant roots seem act as a positive function in root-composite strength, while a negative role in water infiltration. Plant roots have various architectures in different land ecosystems and climatic conditions (Ma et al. 2018). Increasing studies related to soil–root complex focus on the root architectures (Burylo et al. 2011; Li et al. 2016a, b). One major controversy exists, e.g., the plant roots play positive role and enhance slope strength (Arnone et al. 2016). The roots could advance rainfall infiltration, thus contributing an adverse effect on slope stability (Ghestem et al. 2011; Garg et al. 2015). The root–soil composite strength and the hydraulic conductivity are of utmost importance for the rooted soil slope stability.

### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Chapter 2 Analytical Solution to Unsaturated Infiltration**

### **2.1 Introduction**

Rainfall infiltration in unsaturated soil slopes is a classic issue in geotechnical engineering (Conte and Troncone 2012; Iverson 2000; Morbidelli et al. 2018). Factors influencing the soil slope stability due to rainwater infiltration comprise the rainfall characteristics, the saturated permeability coefficient, the geometry of the slope, and the boundary and initial soil moisture conditions (Ali et al. 2014; Wu et al. 2020).

The spatial and temporal evolution of unsaturated infiltration involves a governing partial differential equation that is expressed by the Richards' equation (1931). The equation is highly nonlinear because the hydraulic conductivity and the pressure head depend on matric suction or moisture content. A number of exact and approximated analytical solutions to the infiltration equation were derived in past studies (e.g., Parlange et al. 1997; Basha 2011). Analytical solutions of the linearized infiltration equation were derived as an integral (Chen et al. 2003), as a Laplace transformation (Zhan et al. 2013), and as a Green's function (Basha 1999). While numerical approaches can effectively simulate complex nonlinear infiltrations into an unsaturated soil (Tracy 2006; Wu et al. 2020), analytical solutions can verify these numerical procedures. The incorporation of physically based infiltration expressions better quantifies the infiltration models and causes a more reliable prediction of the water infiltration (Basha 2011). Compared with numerical solutions, the analytical methods are widely used because they can accurately check numerical methods, and concisely represent the pressure head variations with rainfall infiltration (Godt et al. 2012; Wu et al. 2020). Several analytical solutions to Richards' equation were obtained using integral transform methods (Laplace and Fourier algorithms) and others (Qin et al. 2010). The interaction between soil infiltration and displacement was defined as hydro-mechanical coupling effect (Wu et al. 2020). The coupled equation is characterized by strong nonlinearity; thus, linearization is required when solving the equation (Li et al. 2013; Zhu et al. 2022). To derive the analytical solution to infiltration equation considering the coupled hydro-mechanical effect, the soil–water

<sup>©</sup> The Author(s) 2023

L. Wu and J. Zhou, *Rainfall Infiltration in Unsaturated Soil Slope Failure*, SpringerBriefs in Applied Sciences and Technology, https://doi.org/10.1007/978-981-19-9737-2\_2

characteristic curve (SWCC) is represented using an exponential form, and then the analytical solution of pressure head with arbitrary initial conditions can be developed using integral transformation and other methods (Li and Wei. 2018; Wu et al. 2009, 2012, 2016, 2018, 2020).

Many natural slopes are covered with vegetation, which can hydraulically and mechanically affect the slope stability (Lynch 1995; Leung et al. 2017). Vegetation root water uptake changes the pore-water pressure or pressure head of slopes, which is defined as the hydraulic effect of vegetation on slopes (Leung et al. 2015; Ni et al. 2018; Wu et al. 2021). Water is transferred by roots from soil voids to upper stems and leaves. The root water uptake rate and hydraulic conductivities are influenced by a number of factors including the soil type, self-hydraulic resistance, and hydraulic resistance of surrounding soils (Nyambayo and Potts 2010). Based on exponential function SWCC, the analytical solutions for water unsaturated infiltration were developed using Green's function (Ng et al. 2015). The effect of hydrological conditions on the stability of vegetated soil slopes was investigated (Liu et al. 2016; Wu et al. 2022).

This chapter is to derive an analytical solution incorporating both the root contribution and hydro-mechanical coupling. The governing equation considering vegetated root and coupled infiltration and deformation is developed. The analytical solution is obtained using Green's function. Parametric analyses are carried out to investigate the effect of factors on the vegetated slope stability.

### **2.2 1D Analytical Solutions for Unsaturated Seepage**

### *2.2.1 No Coupling*

The one-dimensional infiltration equation can be expressed as (Richards 1931):

$$\frac{\partial}{\partial z} \left[ K(h) \left( \frac{\partial h}{\partial z} + 1 \right) \right] = \frac{\partial \theta}{\partial t} \tag{2.1}$$

where *h* is the pressure head; *z* is the vertical direction; *t* is the rainfall time.

Based on Gardner model (1958) and variable definition (*h*<sup>∗</sup> = e<sup>α</sup>*<sup>h</sup>*−χ in Chap. 1, where α is the desaturation coefficient, χ = e<sup>α</sup>*h*<sup>d</sup> ), the linearized Richards' equation can be obtained:

$$\frac{\partial^2 h^\*}{\partial z^2} + \alpha \frac{\partial h^\*}{\partial z} = c \frac{\partial h^\*}{\partial t} \tag{2.2}$$

where *c* = α(θs − θr)/*K*s.

Tracy (2006) developed an analytical solution to transient infiltration in unsaturated soils:

### 2.2 1D Analytical Solutions for Unsaturated Seepage 17

$$h\_{\mathbf{a}}(z,t) = \frac{1}{\alpha} \ln \left( h\_{\mathbf{l}}^{\*}(z,t) + h\_{\mathbf{s}}^{\*}(z) + \mathbf{e}^{a h\_{\mathbf{d}}} \right) \tag{2.3}$$

in which,

$$h\_s^\*(\mathbf{z}) = \left(1 - \mathbf{e}^{\alpha h\_d}\right) \times \left(1 - \mathbf{e}^{-\alpha \varepsilon}\right) / \left(1 - \mathbf{e}^{-\alpha L}\right) \tag{2.4}$$

$$h\_{\rm t}^{\*}(z,t) = \frac{2\left(1 - e^{a h\_{\rm d}}\right)}{Lc} \mathbf{e}^{a(L-z)/2} \sum\_{m=1}^{\infty} (-1)^{m} \left(\frac{\lambda\_{\rm m}}{\mu\_{\rm m}}\right) \sin(\lambda\_{\rm m} z) \mathbf{e}^{-\mu\_{\rm m} t} \tag{2.5}$$

where λm = *m*π/*L* and μm = ( α<sup>2</sup>/4 + λ<sup>2</sup> m ) /*c*.

For unsaturated slopes (Fig. 2.1), the governing equations are modified as follows:

$$\frac{\partial}{\partial z} \left[ K\_z(h) \left( \frac{\partial h}{\partial z} + \cos \beta \right) \right] = \frac{\partial \theta}{\partial t} \tag{2.6}$$

where β is the slope angle. The analytical solution of the seepage equation along the *z*-axis can be expressed as:

$$h\_{\rm l}^{\*}(z,t) = \frac{2\left(1 - e^{\alpha h\_{\rm d}}\right)}{Lc} \mathbf{e}^{\alpha \cos\beta (L-z)/2} \sum\_{m=1}^{\infty} (-1)^{m} \left(\frac{\lambda\_{\rm m}}{\mu\_{\rm m}}\right) \sin(\lambda\_{\rm m} z) \mathbf{e}^{-\mu\_{\rm m}t} \tag{2.7}$$

$$h\_s^\*(z) = \left(1 - \mathbf{e}^{\alpha h\_d}\right) \frac{1 - \mathbf{e}^{-\alpha \cos \beta z}}{1 - \mathbf{e}^{-\alpha \cos \beta L}}\tag{2.8}$$

**Fig. 2.1** Top and bottom flux boundaries for a soil profile with a finite thickness. The soil layer is between *z* = 0 and *z* = *l*, where *z* is the vertical coordinate: **a** impermeable boundary; and **b**  groundwater level at the bottom

### *2.2.2 Hydro-mechanical Coupling*

Based on Darcy's law, mass and momentum conservation, the equation that governs 1D hydro-mechanical coupling in unsaturated soils, including the consideration of increases in the water table, can be given by (Lloret 1987; Wu et al. 2020):

$$\frac{\partial}{\partial z} \left[ k \frac{\partial}{\partial z} (h + z) \right] = \left( n \beta\_{\text{w}} S\_{\text{r}} \chi\_{\text{w}} + n \frac{\partial S\_{\text{r}}}{\partial h} \right) \frac{\partial h}{\partial t} - S\_{\text{t}} \alpha\_{\text{c}} \frac{\partial \varepsilon\_{\text{z}}}{\partial t} \tag{2.9}$$

$$\frac{\partial}{\partial z} \left[ E\varepsilon\_{\rm v} - \frac{E}{F} (\mu\_{\rm a} - \mu\_{\rm w}) \right] + [nS\_{\rm r}\rho\_{\rm w} + (1 - n)\rho\_{\rm s}]\mathbf{g} = 0 \tag{2.10}$$

where εv (εv = ε*z* for one-dimensional problem) is the total volumetric strain of the soil mass, which is greater than zero during compression and less than zero during swelling; *n* is the percentage of voids; *S*r is the saturation; (*u*a − *u*w) is the matric suction; αc is the hydro-mechanical coupling coefficient (0 ≤ αc ≤ 1), which is determined by the bulk modulus of the solid skeleton and the bulk modulus of the solid soil (Wu et al. 2020); *E* is the elastic modulus of the soil that takes account of alterations in the net normal stresses (Wu et al. 2020); *F* is the elastic modulus of the soil due to variations in matric suction, which is assumed to be a function of stress; βw is the compressibility of the fluid; ρw is the density of water; ρs is the density of the soil phase; and *g* is the acceleration of gravity.

Equations (2.9) and (2.10) govern the coupling of 1D deformation and seepage in unsaturated soils and account for the groundwater table changes.

On the assumption that pore-air pressure in the soil mass is kept constant; the derivative of Eq. (2.10) with respect to *t* can be written as:

$$\frac{\partial \varepsilon\_{\rm V}}{\partial t} = -\frac{\varkappa\_{\rm w}}{F} \frac{\partial h}{\partial t} \tag{2.11}$$

From Eqs. (2.9) and (2.11), the following equation can be derived:

$$\frac{\partial}{\partial z} \left[ k \frac{\partial}{\partial z} (h + z) \right] = \left( \beta\_{\text{w}} \frac{\partial \rho}{\partial h} + n \frac{\text{d}S\_{\text{r}}}{\text{d}h} \right) \frac{\partial h}{\partial t} + \frac{\wp\_{\text{w}} S\_{\text{r}} \alpha\_{\text{c}}}{F} \frac{\partial h}{\partial t} \tag{2.12}$$

Based on assumption that the coefficient of permeability and the moisture content vary exponentially with the pore-water pressure head, that is, *k*(*h*) = *k*se<sup>α</sup>*<sup>h</sup>*and θ(*h*) = θse<sup>α</sup>*<sup>h</sup>*(*h* < 0, *t* > 0). Here, *k*s is the saturated hydraulic conductivity, θs is the volumetric moisture at saturation. Because *S*r = θ(*h*)/θs when the fluid compressibility is neglected (βw = 0), Eq. (2.12) can be written as follows:

$$\frac{\partial}{\partial z} \left[ k \frac{\partial}{\partial z} (h+z) \right] = M \mathbf{e}^{ah} \frac{\partial h}{\partial t} \quad (t>0) \tag{2.13}$$

where *M* = θsα + γ <sup>w</sup>αc/*F*.

With dimensionless variables *Z* = α*z* and *T* = α2*k*s*t*/*M* introduced into a function of *W*(*Z*, *T*) = e<sup>α</sup>*<sup>h</sup>*· e*<sup>Z</sup>*/2+*<sup>T</sup>*/4, Eq. (2.13) can be rewritten as:

$$\frac{\partial W}{\partial T} = \frac{\partial^2 W}{\partial Z^2} \tag{2.14}$$

The boundaries comprise a base and top one in Fig. 2.1. In the literature of analytical solutions, the base boundary was usually assumed to coincide with a stationary groundwater table and the pressure head was set zero (Wu et al. 2020). However, in this book, a zero flux is considered at the base boundary. The hydraulic boundary condition is given by:

$$h\Big|\_{z=0} = 0\tag{2.15}$$

or

$$k\frac{\partial h}{\partial z} + k\bigg|\_{z=0} = 0\tag{2.16}$$

The top boundary in Fig. 2.1 is controlled by pressure head or rainfall intensity (*q*) at the ground surface, and it is written as:

$$h\Big|\_{z=l} = h\_0\tag{2.17}$$

or

$$k\frac{\partial h}{\partial z} + k\Bigg|\_{z=l} = q(t) \tag{2.18}$$

in which, *l* is the depth in the 1D unsaturated infiltration model.

Based on a Fourier integral transformation (Ozisik 1989), the exact solution to Eq. (2.14) can be derived considering different boundaries (Wu et al. 2020).

### **2.3 2D Analytical Solutions of Rainfall Infiltration in Unsaturated Soils**

The 2D Richards' equation in mixed format is expressed as:

$$\frac{\partial}{\partial x}\left[K\_x(h)\frac{\partial h}{\partial x}\right] + \frac{\partial}{\partial z}\left[K\_z(h)\left(\frac{\partial h}{\partial z} + 1\right)\right] = \frac{\partial \theta}{\partial t} \tag{2.19}$$

**Fig. 2.2** Schematic diagram of 2D transient seepage model in unsaturated soils

in which, *Kx* and *Kz* are the hydraulic conductivities along *x*- and *z*-directions, respectively.

Similarly, the 2D linearized Richards' equation can be written as:

$$\frac{\partial^2 h^\*}{\partial x^2} + \frac{\partial^2 h^\*}{\partial z^2} + \alpha \frac{\partial h^\*}{\partial z} = c \frac{\partial h^\*}{\partial t} \tag{2.20}$$

The mathematical model is shown in Fig. 2.2. *L* and *W* represent the height and length, respectively. The normalized boundary conditions are as follows:

$$h(0, z, t) = h\_{\mathbb{d}} \tag{2.21}$$

$$h(W, z, t) = h\_{\mathbb{d}} \tag{2.22}$$

$$h(\mathbf{x}, \mathbf{0}, t) = h\_{\mathbf{d}} \tag{2.23}$$

$$h(\mathbf{x}, L, t) = \ln\left( \left( 1 - \mathbf{e}^{ah\_d} \right) \sin\left(\frac{\pi x}{W}\right) + \mathbf{e}^{ah\_d} \right) / \alpha \tag{2.24}$$

The normalized analytical solutions *h*<sup>∗</sup> <sup>t</sup>(*x*,*z*, *t*) and *h*<sup>∗</sup> <sup>s</sup>(*x*,*z*, *t*) for this twodimensional model and can be expressed as follows (Tracy 2006):

$$h\_{\rm t}^{\*}(\mathbf{x}, z, t) = \frac{2\left(1 - \mathbf{e}^{ah\_d}\right)}{Lc} \sin\left(\frac{\pi x}{W}\right) \mathbf{e}^{a(L-z)/2} \sum\_{m=1}^{\infty} (-1)^m \left(\frac{\lambda\_{\rm m}}{\mu\_{\rm m}}\right) \sin(\lambda\_{\rm m} z) \mathbf{e}^{-\mu\_{\rm m} t} \tag{2.25}$$

$$h\_s^\*(x, z, t) = \left(1 - \mathbf{e}^{a h\_4}\right) \mathbf{e}^{\frac{a(l-\varepsilon)}{2}} \sin\left(\frac{\pi x}{W}\right) \frac{\sinh(\beta\_1 z)}{\sinh(\beta\_1 L)}\tag{2.26}$$

where β1 <sup>=</sup> <sup>√</sup> α<sup>2</sup>/4 + (π/*W* ) 2 .

### **2.4 Analytical Solution of Water Infiltration of Vegetated Slope Considering the Coupling Effects**

Coupling between water infiltration and mechanical deformation in unsaturated soils is central to many natural and man-made systems in civil and environmental engineering. During water infiltration, the pore-water pressure or pressure head is redistributed, on one hand by the hydraulic properties of the unsaturated soils including retention characteristics and permeability, and, on the other hand, by the external loading due to climate conditions (rainfall intensity, duration, and evapotranspiration rate). Changes in the pore-water pressure or pressure head are generated by infiltration, which in turn modifies the hydraulic domain and induces deformations of the unsaturated soils. Alternatively, any variation in the mechanical loading can exert an effect on the infiltration process. It is indeed the hydro-mechanical coupled response of an unsaturated soil that is responsible for the most common instabilities associated with water infiltration: landslides and settlements, due to collapse or shear strength reduction (Thorel et al. 2011; Wu et al. 2020).

Recently, ecological protection technologies have become popular for slope environmental restoration and treatment (Tan et al. 2019; Broda et al. 2020). The stability analysis of vegetated slopes is a hot point in geotechnical engineering. Several studies have analyzed the stability of infinite vegetated slopes (Feng et al. 2020). However, few studies have been reported on the analytical solution of vegetated slope stability considering the hydro-mechanical coupling. Developing analytical solutions for water infiltration in unsaturated soil slopes is a significant issue in practical engineering. The main objective of this section is to derive an analytical solution considering both the root effect and hydro-mechanical behavior. The governing equation considering vegetated root and coupled infiltration-deformation is derived. The analytical solution is developed using Green's function method, which is then compared with the numerical solution. Parametric analyses are performed to investigate the effect of factors on the infinite vegetated slope stability.

### *2.4.1 Governing Equations of Rooted Unsaturated Soils*

To simplify the problem, the following assumptions are made:


An unsaturated soil slops with vegetation is shown in Fig. 2.3, *L*\* is the thickness of a soil (m), *L*<sup>∗</sup> 1 is the depth of the rooted zone (m), *L*<sup>∗</sup> 2 is the thickness of the unrooted zone (m). The Richards' equation in a reference coordinate system (*oxz*) can be expressed as:

$$\frac{\partial}{\partial x}\left[k(h)\frac{\partial h}{\partial x}\right] + \frac{\partial}{\partial z}\left[k(h)\frac{\partial h}{\partial z} + k(h)\right] = C(h)\frac{\partial h}{\partial t} \tag{2.27}$$

where *h* is the pressure head (m); *k*(*h*) is the hydraulic conductivity (m/s); *C*(*h*) = dθ/d*h* is the differential water capacity (m−1); θ is the volumetric water content; *t*  is the infiltration duration (s). Equation (2.27) needs to be modified to be suitable for rooted soil slopes. The transformation relationship between the slope coordinate system (*ox*\**z*\*) and the reference coordinate system (*oxz*) is described as follows:

$$\mathbf{x}^\* = \mathbf{x}\cos\beta - z\sin\beta \tag{2.28a}$$

$$z^\* = x\sin\beta + z\cos\beta\tag{2.28b}$$

**Fig. 2.3** Graphical representation of rainfall-induced landslides

Substituting Eq. (2.28) into Eq. (2.27), according to Assumption (1), the improved Richards' equation for rainfall infiltration into soil slopes can be given by:

$$\frac{\partial}{\partial z^\*} \left[ k(h) \frac{\partial h}{\partial z^\*} \right] + \frac{\partial k(h)}{\partial z^\*} \cos \beta = C(h) \frac{\partial h}{\partial t} \tag{2.29}$$

The sink term proposed by Raats (1979) is used to describe the water uptake of roots:

$$S(z^\*) = g(z^\*)T\tag{2.30}$$

where *g*(*z*\*) is the shape function of roots at depth *z*\* ([m/s]−1); and *T* is the transpiration rate (m/s), which is affected by weather and leaf area index. Substituting Eq. (2.30) into Eq. (2.29), according to Assumption (2), the modified Richards equation considering the water uptake by vegetation roots can be obtained (Wu et al. 2022):

$$\frac{\partial}{\partial z^\*} \left[ k(h) \frac{\partial \psi}{\partial z^\*} \right] + \frac{\partial k(h)}{\partial z^\*} \cos \beta - S(z^\*) [z^\* - L\_2^\*] = C(h) \frac{\partial h}{\partial t} \tag{2.31}$$

where ⧼ *z*<sup>∗</sup> − *L*<sup>∗</sup> 2 ⧽ = ( *z*<sup>∗</sup> − *L*<sup>∗</sup> <sup>2</sup>, *z*<sup>∗</sup> ≥ *L*<sup>∗</sup> 2 0, *z*<sup>∗</sup> < *L*<sup>∗</sup> 2 is the sign function.

Equation (2.31) ignores the hydraulic coupling in unsaturated soils. According to the law of conservation of mass, the strict expression of the term *C*(*h*) <sup>∂</sup>*<sup>h</sup>* <sup>∂</sup>*<sup>t</sup>*on the right side of Eq. (2.31) is <sup>1</sup> ρ ∂ <sup>∂</sup>*<sup>t</sup>*(ρ*nS*r) (Kim 2000; Wu et al. 2020, 2022), where ρ, *n*, and *S*<sup>r</sup> represent the density of water, soil porosity, and degree of saturation. *S*r = (θ(*h*) − θr)/(θs − θr), θs is the saturated water content, and θr is the residual water content. Then one can obtain:

$$\frac{\partial}{\partial t}(\rho n S\_{\rm w}) = \theta \frac{\partial \rho}{\partial t} + \rho \frac{\partial S\_{\rm r}}{\partial t} n + \rho S\_{\rm r} \frac{\partial n}{\partial t} \tag{2.32}$$

the first term on the right side of Eq. (2.32) represents the water compression term, which is 0 according to Assumption (3); the second term is equivalent to *C*(*h*) <sup>∂</sup>*<sup>h</sup>* ∂*t* ; the third term represents the soil skeleton term, and ρ*S*<sup>r</sup> ∂*n*  <sup>∂</sup>*<sup>t</sup>*= ρ*S*rα<sup>c</sup> ∂ε<sup>v</sup> <sup>∂</sup>*<sup>t</sup>*, where αc is Biot's hydro-mechanical coupling parameter (0 < αc ≤ 1), Ev is volumetric strain. Therefore, Eq. (2.31) can be rewritten as:

$$\frac{\partial}{\partial z^\*} \left[ k(h) \frac{\partial h}{\partial z^\*} \right] + \frac{\partial k(h)}{\partial z^\*} \cos \beta - S(z^\*) \langle z^\* - L\_2^\* \rangle = C(h) \frac{\partial h}{\partial t} - S\_t a\_\varepsilon \frac{\partial \varepsilon\_\mathbb{V}}{\partial t} \quad (2.33)$$

Substituting Eq. (2.11) into Eq. (2.33), the governing equation considering the hydro-mechanical coupling in unsaturated soil slopes can be obtained:

24 2 Analytical Solution to Unsaturated Infiltration

$$\frac{\partial}{\partial z^\*} \left[ k(h) \frac{\partial h}{\partial z^\*} \right] + \frac{\partial k(h)}{\partial z^\*} \cos \beta - S(z^\*) [z^\* - L\_2^\*] = C(h) \frac{\partial h}{\partial t} + \frac{\chi\_\text{w} (\theta - \theta\_\text{r}) a\_\text{c}}{(\theta\_\text{s} - \theta\_\text{t}) F} \frac{\partial h}{\partial t} \tag{2.34}$$

### *2.4.2 Analytical Solutions*

According to Assumption (1), the bottom and surface boundary conditions of the slope are written as, respectively:

$$h(0, t) = 0 \quad t > 0\tag{2.35}$$

$$\left[k(h)\frac{\partial h(L^\*,t)}{\partial z^\*} + k(h)\cos\beta\right] = q\cos\beta \quad t>0\tag{2.36}$$

where *q* is the rain intensity (m/s).

The hydraulic conductivity and the volumetric water content of unsaturated soils can be expressed as (Gardner 1958):

$$k(h) = k\_s \mathbf{e}^{\alpha h} \tag{2.37}$$

$$\theta(h) = \theta\_{\mathbf{r}} + (\theta\_{\mathbf{s}} - \theta\_{\mathbf{r}}) \mathbf{e}^{ah} \tag{2.38}$$

where *k*s represents the saturated hydraulic conductivity (m/s); α is the desaturation coefficient (kPa−1). Substituting Eq. (2.38) into Eq. (2.34), one has:

$$\frac{\partial}{\partial z^\*} \left[ k \frac{\partial h}{\partial z^\*} \right] + \frac{\partial k}{\partial z^\*} \cos \beta - S(z^\*) \langle z^\* - L\_2^\* \rangle = M \mathbf{e}^{ah} \frac{\partial h}{\partial t} \tag{2.39}$$

where *M* = (θs − θr)α + γ <sup>w</sup>αc/*F*.

The soil–water uptake function (Eq. 2.30) can be simplified as (Lynch 1995):

$$S(z^\*) = T/L\_1^\* \tag{2.40}$$

Substituting Eqs. (2.37), (2.38), and (2.40) into Eq. (2.39) leads to:

$$\frac{\partial^2 k}{\partial z^{\*2}} + \frac{\partial k}{\partial z^\*} \alpha \cos \beta - T / L\_1^\* [z^\* - L\_2^\*] = \frac{M}{K\_s} \frac{\partial k}{\partial t} \tag{2.41}$$

Here, the variables are defined as:

$$\begin{cases} Z = z^\* \cos \beta \\ H = L^\* \cos \beta \\ H\_1 = L\_1^\* \cos \beta \\ H\_2 = L\_2^\* \cos \beta \\ K = k/k\_s \\ Q\_\sf a = q\_\sf a / k\_s \\ Q\_\sf b = q\_\sf b / k\_s \\ P = k\_s/M \end{cases} \tag{2.42}$$

where *q*a and *q*b represent the previous rain intensity (m/s) and current rain intensity (m/s).

Substituting Eq. (2.42) into Eq. (2.41) leads to:

$$\frac{\partial^2 K}{\partial Z^2} + \frac{\partial K}{\partial Z} \alpha - \frac{S(z^\*/\cos \beta)[z^\*/\cos \beta - L\_2^\*/\cos \beta]}{k\_s \cos^2 \beta} = \frac{1}{P \cos^2 \beta} \frac{\partial K}{\partial t} \tag{2.43}$$

By combining Eqs. (2.35), (2.36), and (2.43), the analytical solution of pressure head can be obtained using Green's function as follows:

$$K\_{\rm stc} = \begin{cases} \exp(-\alpha Z) + Q\_{\rm a} [\exp(-\alpha Z) - 1] \\ + \frac{T}{K\_s \cos \beta H\_1} [\exp(-\alpha Z) - 1](H - H\_1), & Z < H\_2 \\ \exp(-\alpha Z) + Q\_{\rm a} [\exp(-\alpha Z) - 1] + \frac{T}{K\_s \cos \beta H\_1} \\ \left\{ \begin{bmatrix} \exp(-\alpha Z) - 1 \\ Z - H\_1 - \alpha^{-1} \exp(\alpha Z) + \alpha^{-1} \exp(\alpha H\_1) \end{bmatrix} \right\}, & Z \ge H\_2 \end{cases} \tag{Stady-state}$$

$$\begin{split} K &= K\_{\text{ste}} + 8 \frac{\alpha P}{k\_{\text{s}}} \cos^{2} \beta \exp\left[\frac{\alpha (H - Z)}{2}\right] \\ &\times \sum\_{i=1}^{\infty} \frac{\left[\lambda\_{i}^{2} + 0.25 \alpha^{2}\right] \sin(\lambda\_{i} H) \sin(\lambda\_{i} Z)}{2\alpha + \alpha^{2} H + 4H \lambda\_{i}^{2}} G(t) \qquad (\text{Transient}) \end{split} \tag{2.45}$$

where

$$\mathcal{G}(t) = \int\_0^t (\mathcal{Q}\_\mathfrak{a} - \mathcal{Q}\_\mathfrak{b}) k\_\text{s} \exp[-P \cos^2 \beta \left(\lambda\_i^2 + 0.25 \alpha^2\right)(t - \tau)] d\tau \tag{2.46}$$

λ*i* is the *i*th root of the transcendental equation tan(λ*H*) + (2λ/*a*) = 0. The transient solution changes with time. The "steady-state" infiltration evolves over formally infinite time, and thus the "steady-state" solution isn't related to time. However, transient infiltration depends on time, and the analytical solution to transient infiltration is called transient solution. Rainfall infiltration into soil slopes is commonly a transient issue, thus transient solutions are greatly meaningful.

When *Q*b is a constant, Eq. (2.45) can be simplified as:

$$\begin{split} K &= K\_{\text{ste1}} - 8(Q\_{\text{a}} - Q\_{\text{b}})\alpha \cos^{2}\beta \exp\left[\frac{\alpha(H - Z)}{2}\right] \\ &\times \sum\_{i=1}^{\infty} \frac{\sin(\lambda\_{i}H)\sin(\lambda\_{i}Z)\exp\left[-P\cos^{2}\beta \left(\lambda\_{i}^{2} + 0.25\alpha^{2}\right)t\right]}{2\alpha + \alpha^{2}H + 4H\lambda\_{i}^{2}} \end{split} \tag{2.47}$$

where

$$K\_{\rm site} = \begin{cases} \exp(-\alpha Z) + Q\_b[\exp(-\alpha Z) - 1] \\ + \frac{T}{k\_s \cos \beta H\_1}[\exp(-\alpha Z) - 1](H - H\_1), & Z < H\_2 \\ \exp(-\alpha Z) + Q\_b[\exp(-\alpha Z) - 1] + \frac{T}{k\_s \cos \beta H\_1} \\ \left\{ \begin{bmatrix} \exp(-\alpha Z) - 1 \\ Z - H\_1 - \alpha^{-1} \exp(\alpha Z) + \alpha^{-1} \exp(\alpha H\_1) \end{bmatrix} \right\}, & Z \ge H\_2 \end{cases} (2.48)$$

The pore-water pressure head is calculated as follows:

$$h = \frac{\ln(K)}{\alpha} \tag{2.49}$$

### *2.4.3 Comparison of Analytical and Numerical Solutions*

The finite element method is employed to obtain the numerical solutions to Eqs. (2.34)–(2.36). In this book, COMSOL software is employed to implement numerical solutions, which has been used in geotechnical, hydraulic, and civil engineering.

The parameters adopted here are listed in Table 2.1, which includes the hydraulic parameters (*k*s, α, θs, and θr) (Liu et al. 2016). The parameters describe the hydraulic properties of loess, and the parameters of vegetation roots are determined according to the statistics of shrubs (Feng et al. 2020).

The analytical and numerical solutions of pressure head at *t* = 0 h, 20 h, and 30 h are shown in Fig. 2.4. The analytical solutions at *t* = 0 h are obtained by Eq. (2.44), and those at *t* = 20 h and 30 h are obtained by Eq. (2.45). According to the root mean square error RMSE and coefficient of determination (*R*2) shown in Fig. 2.4, it can be seen that RMSE between the numerical and analytical solutions is less than 0.005 m, and *R*2 is very close to 1. That is, the error between the numerical and analytical solutions is very small. Compared with the numerical solution, the analytical solution concisely describes the pressure head variations with rainfall infiltration.


# *2.4.4 Parametric Analyses*

Factor of safety (*F*s) is an important indicator used to assess the stability of slopes, and slopes with *F*s greater than 1 are generally considered to be stable, and vice versa. *F*s is calculated based on the limit equilibrium method and strength theory of rooted unsaturated soils (Fredlund and Rahardjo 1993; Ku et al. 2018), as follows:

$$F\_s = \frac{c + c\_\mathbf{r}}{(L^\* - z^\*)\chi \cos \beta \sin \beta} + \frac{\tan \varphi}{\tan \beta} - \frac{\chi\_\mathbf{w} h[(\theta - \theta\_\mathbf{r})/(\theta\_\mathbf{s} - \theta\_\mathbf{r})] \tan \varphi}{(L^\* - z^\*)\chi \cos \beta \sin \beta} \tag{2.50}$$

where *c*r is the root cohesion (kPa), obtained by Wu–Waldron criterion (Wu 1976), i.e., *c*r = ς*T*s, *T*s is the root tensile strength (kPa), ς is the ratio of root cross-sectional area, which is 0.00025 in this book (Leung 2014). *c* and ϕ are the soil effective cohesion and effective frictional angle, while the rooted soil strength contributes to root cohesion. In this section, the effects of slope angle, rainfall intensity, transpiration rate, and suction-based modulus of elasticity on *F*s and *F*s ratio will be investigated. The *F*s ratio, defined as the ratio of *F*s during rainfall to the initial value (*t* = 0), shows the variations in *F*s during rainfall and is also an important parameter for analyzing rainfall-include landslides.

### **2.4.4.1 Effect of Slope Angle**

Here, the effects of slope angle on the pressure head, *F*s and *F*s ratio (defined as the ratio of *F*s during rainfall to the initial value (*t* = 0)) are investigated. The parameters are as follows: Rain intensity is 0.7*k*s, the tensile strength of the roots is 10 MPa (this parameter is determined by the statistics of shrubs (Leung 2014)), the unit weight of the soil (γ ) is 20 kN/m3, the soil effective cohesion and effective frictional angle are 18 kPa and 28°, and other parameters are listed in Table 2.1.

Figure 2.5 represents the pressure head in slopes with slope angles of 10, 30, and 50° during rainfall. The conclusions can be drawn as follows (Fig. 2.5): (i) the pressure head of slopes decreases during rainfall; (ii) the pressure head in the shallow slope is more sensitive than that in the deep slope during rainfall; (iii) the smaller the slope angle, the larger the pressure head at the same position, which is mainly caused by the initial pressure head profile (the initial pressure head of the slope is negatively correlated with slope angle); (iv) the smaller the slope angle, the greater the variation in the pressure head at the same time, this is because the small slope angle has a positive effect on rainfall infiltration.

*F*s (related to depth) with slope angles of 10, 30, and 50° during rainfall are described in Fig. 2.6a. Figure 2.6a states the following points: (a) *F*s of the shallow slope is greater than that of the deep slope; (b) *F*s decreases during rainfall, which can be explained in Fig. 2.5; (c) increase at *z* = 4.5 m due to the reinforcement effect of vegetation roots. *F*s ratio *F*s ratio with slope angles of 10, 30, and 50° during rainfall are represented in Fig. 2.6b. The larger the slope angle, the larger *F*s ratio in Fig. 2.6b. This demonstrates that the larger the slope angle, the smaller *F*s. In Fig. 2.6a, *F*s suddenly the stability of a gentle slope is more sensitive to rainfall than that of the deep slope.

### **2.4.4.2 Effect of Rain Intensity**

The variations in the dimensionless rainfall intensity over the duration of the rainfall event for both coupled and uncoupled analyses are shown in Fig. 2.7. The parameters are *k*s = 10−5 m/s, θs = 0.4, |*F*| = 5 × 103 kPa, and α = 0.01 cm−1 (Van Genuchten 1980). The model height was 400 cm (Fig. 2.7). The value of *F* can be positive or negative. A negative *F* means an expansive soil where a suction decrease leads to soil volume increase, while a positive *F* denotes a collapsible soil where a soil suction decrease leads to soil volume decrease (Wu et al. 2020). Compared with the bottom boundary of a stationary groundwater table, the impermeable bottom boundary leads to more pronounced coupling effects in the lower part of the soil layer. It is also noted that the groundwater ponding occurs at the bottom boundary (*t* = 15 h, Fig. 2.7), particularly for an expansive soil (*F* < 0). Waterfall intensity plays a significant role in the advancement of the wetting front, and the pressure head profile moves more quickly as the rainfall intensity increases. The coupling effect is also closely linked with the rainfall intensity (Wu et al. 2020).

Here, the effect of rain intensity on the pressure head, *F*s and *F*s ratio are investigated. Figure 2.8 describes the pressure head of the slope with rainfall intensities of 0.5*k*s, 0.7*k*s and 0.9*k*s. The remaining parameters are listed in Table 2.1. In Fig. 2.8, the smaller the rainfall intensity, the larger the pressure head at the same position.

Figure 2.9a, b describe *F*s, and *F*s ratio of the slope with rain intensities of 0.5*k*s, 0.7*k*s, and 0.9*k*s, respectively. The variation in *F*s with time and depth is similar to that in Fig. 2.6a. The larger the rain intensity, the smaller *F*s or *F*s ratio of the slope. The stability of the shallow slope is more sensitive to that the deep slope under conductivity is an important parameter of soil seepage capacity, different rainfall intensities (Fig. 2.6). Saturated hydraulic conductivity is an important parameter of soil seepage capacity. Figure 2.10 represents the influence of saturated hydraulic

conductivity on *F*s and *F*s ratio of rooted soil slopes. With increasing saturated hydraulic conductivity, *F*s and *F*s ratio at the same depth decrease.

### **2.4.4.3 Effect of Transpiration Rate**

The transpiration rate is an important parameter. The effect of transpiration rate on the pressure head, *F*s and *F*s ratio, is investigated here. The parameters are as follows: Rain intensity is 0.7*k*s, and other parameters are the same as those in Sect. 2.4.4.2.

Figure 2.11 represents the pressure head of unsaturated soil slopes with transpiration rates of 3, 4.5, and 6 mm/d. Transpiration rates depend on the vegetation and environmental conditions, and the transpiration rates are determined according to the statistics of shrubs (Leung and Ng 2013). The pressure head varies over time, and depth is the same as that in Sect. 2.4.4.1. The smaller the transpiration rate, the smaller the pressure head at the same position. The transpiration of vegetation reduces water content in slopes, causing an increase in soil suction. Figure 2.12a depicts *F*<sup>s</sup> of slopes with transpiration rates of 3, 4.5, and 6 mm/d, respectively. The larger the transpiration rate, the larger safety factor of soil slopes. This is because the vegetation root uptake water reduces the water content of soil slopes. Figure 2.12b represents *F*s ratio of slopes with transpiration rates of 3, 4.5, and 6 mm/d, respectively. The larger the transpiration rate, the smaller *F*s ratio at the same position.

### **2.4.4.4 Effect of the Suction-Based Modulus of Elasticity**

The suction-based modulus of elasticity (*F*) is key for hydro-mechanical coupling. The effect of suction-based modulus of elasticity on the pressure head, *F*s, and *F*s

ratio, was investigated here. Three cases (*F* = 103 kPa, − 102 kPa, and no hydromechanical coupling) are considered here. The governing equation without the hydromechanical coupling is Eq. (2.31), and the boundary conditions are the same as those considering the coupling effect. The transpiration rate in this section is 4.5 mm/d, and the other parameters are kept unchanged as Sect. 2.4.4.3.

Figure 2.13 represents the pressure head in soils with different suction-based moduli of elasticity. The variations of pressure head with suction-based modulus of elasticity are summarized as follows: (i) When *F* is positive, the hydro-mechanical coupling causes the increase of pore-water pressure head in the slope. Water flow in expansive soils is faster than that in collapsible soils. (ii) The effect of hydromechanical coupling on the pressure head in the slope becomes obvious with the decrease of the absolute value *F*.

Figure 2.14a represents *F*s of soil slopes with different suction-based moduli of elasticity. Variation in safety factor with *F* is as follows: (i) When *F* is positive (negative), the hydro-mechanical coupling effect results in increases (decreases) of *F*s. (ii) The effect of hydro-mechanical coupling on *F*s becomes marked with decreasing absolute value of *F*.

Figure 2.14b describes *F*s ratios of soil slopes with different suction-based moduli of elasticity. Figure 2.14b clearly demonstrates that small absolute value of *F* will significantly reduce the factor of safety of slopes.

### **2.5 Discussions and Conclusions**

### *2.5.1 Discussions*

*F*s ratios caused by the hydro-mechanical coupling effect of bare slopes (transpiration rate is 0) and rooted slopes (transpiration rate is 4.5 mm/d) are compared. The

evaporation of bare soils caused by temperature is ignored. The *F*s ratios of the bare and rooted slopes at *t* = 20 h and 40 h are described in Fig. 2.15. The parameters are the same as those in Sect. 2.4.4.4. In Fig. 2.15a, *F*s ratios of the rooted slope are smaller than those of the bare slope. The stability changes caused by the hydromechanical coupling in rooted soil slopes are greater than those of bare slopes. A similar conclusion can be drawn from Fig. 2.15b.

It should be emphasized that the proposed analytical solution has two main limitations: (i) This solution is only suitable for rooted slopes with uniform root architecture. Existing results (e.g., Ng et al. 2015; Liang et al. 2017) have indicated that vegetation roots have complex shape, including exponential root, triangular root, and uniform root. The book tries to obtain analytical solutions considering different root architecture functions in future work. (ii) The boundary condition of soil slopes is transformed from flow boundary to pressure head boundary due to saturated hydraulic conductivity less than rain intensity, which may affect the results of analytical solutions. Although there still occur some limitations, the proposed analytical solution incorporates the root effect and hydro-mechanical coupling for the first time. In actual engineering, the proposed analytical solution can be easily applied to examine rainfall-induced landslides in rooted soil regions.

**Fig. 2.15** *F*s ratios of the bare slope and vegetated slope: **a** *t* = 20 h, **b** *t* = 40 h

Figure 2.16 represents the variations in the dimensionless rainfall intensity over the duration of the rainfall event for both the coupled and uncoupled conditions. The parameters used are *k*s = 10−5 m/s, θs = 0.4, |*F*| = 5 × 103 kPa, and α = 0.01 cm−<sup>1</sup> (Van Genuchten 1980; Wu et al. 2020). The value of *F* can be positive or negative; a negative *F* means an expansive soil where a suction decrease leads to soil volume increase, while a positive *F* denotes a collapsible soil where a suction decrease leads to soil volume decrease (Wu et al. 2020). When *t* = 1 h, the wetting front moves to a depth of 120 cm in both the uncoupled analysis and coupled analysis with *F* > 0. However, with *F* < 0 the wetting front reaches 180 cm in the coupled analysis. When *t* = 1 h, the difference in the pressure head for the coupled (*F* < 0) and uncoupled conditions is 27.5 cm at the upper boundary, and the maximum difference in the pressure head within the unsaturated zone is 75.1 cm at a soil depth of 340 cm. When *t* = 15 h, the difference in the pressure head at the top boundary between the coupled (*F* > 0) and uncoupled conditions is 13.7 cm, and the maximum difference in the pressure head within the soil is 61.8 cm at a depth of 160 cm. The coupling effect becomes more apparent with increasing time, particularly in the bottom part

**Fig. 2.16** Pore-water pressure profile for different boundaries

of the soil layer. Compared with the case where the base boundary coincides with the stationary groundwater table (Wu et al. 2020), the case with zero flux at the base boundary leads to more pronounced coupling effects in the lower part of the soil column. It is also noted that groundwater ponding occurs at the base boundary (*t* = 15 h, Fig. 2.16).

As the infiltration time increases and the wetting front moves downward, the pressure head increases rapidly in the shallow depths of the unsaturated soils. Rainfall intensity plays an important role in the advancement of the wetting front, and the pressure head profile moves more quickly as the rainfall intensity increases. The coupling effect is also closely linked with the rainfall intensity. Other studies reported that the coupling effect is more noticeable in a shallow layer of an unsaturated soil (Wu et al. 2020). However, in this book, the coupling effect becomes more noticeable in the base layer, particularly with increasing time. As the rainfall accumulates at the base boundary, the coupled effect becomes more apparent there. Figure 2.16 describes the effect of boundary on the pressure profile considering the coupling effect. The boundary effect is more marked than the coupling. The zero-flux base boundary also leads to groundwater ponding in the lower part of the soils. The pressure head profile moves more quickly as the thickness of the soil layer decreases.

### *2.5.2 Conclusions*

The governing equation considering hydro-mechanical coupling of unsaturated rooted slopes is derived. The analytical solution is obtained using Green's function method. This solution is compared with the finite element method. The parametric studies were carried out to investigate the effect of slope angle, rain intensity, suctionbased elastic modulus, and transpiration rate on the pore-water pressure head and slope safety factor. The following conclusions can be obtained:


### **References**


Wu LZ, Xu Q, Wang T (2018) Incorporating hydromechanical coupling of un-saturated soils into the analysis of rainwater-induced groundwater ponding. Int J Geomech 18(6): 06018010

Wu LZ, Huang RQ, Li X (2020) Hydro-mechanical analysis of rainfall-induced landslides. Springer


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Chapter 3 Numerical Solutions to Infiltration Equation**

### **3.1 Introduction**

Unsaturated infiltration issues occur in many fields, such as rainfall-induced soil slope failures (Wu et al. 2020a, b; Jiang et al. 2022), solute migration simulation (Cross et al. 2020), and coal seam water injection and coalbed methane extraction (Wang et al. 2020). Among them, the Richards' equation (Richards 1931) is the basic governing equation for the numerical simulation of unsaturated infiltration. The effective and reliable numerical solution of the Richards' equation is of great significance to scientific research and production in related fields. Before the advent of computers, the investigations on unsaturated infiltration mainly focused on the analytical solution of the infiltration equation, that is, the method of directly solving differential equations using relevant mathematical means. The studies on analytical solutions under certain conditions still have very important theoretical and practical significance. Analytical solutions to the Richards' equation can be obtained under some simplified conditions (Broadbridge et al. 2017). For example, Srivastava and Yeh (1991) employed an exponential model describing the soil-water character curve (SWCC) to derive transient analytical solutions for one-dimensional infiltration in homogeneous and layered soils. Parlange et al. (1999) solved a complex series solution of the one-dimensional Richards equation in order to obtain the infiltration flux. Tracy (2006) proposed two-dimensional and three-dimensional analytical solutions to rainfall infiltration into homogeneous soils based on exponential function. Wu et al. (2016, 2020a) used the Laplace transform method to obtain an analytical solution considering the coupling of infiltration and deformation in unsaturated soils during rainfall, and used it to analyze the stability of infinite unsaturated slopes due to rainfall infiltration. However, the permeability coefficients and soil–water characteristic curves are very complex in reality, which leads to highly nonlinear partial differential equations analytical solutions of which are very difficult to obtain. Therefore, numerical methods are often developed to solve the Richards' equation under

common conditions (Zha et al. 2017; Ku et al. 2018; Zeng et al. 2018; Zhu et al. 2019; Eini et al. 2020).

With the development of computer techniques, numerical methods, namely finite difference method (FDM) (Liu et al. 2017), finite element method (FEM) (Crevoisier et al. 2009), and the finite volume method (FVM) (Liu 2017; Eymard et al. 2006), are increasingly used in infiltration analysis. The discretization of the time derivative is usually performed using the backward difference method (Pop and Schweizer 2011). For example, Patankar (1980) summarized FVM for numerical dis-cretization of heat transfer and fluid flow. Wang and Anderson (1982) introduced the application of finite difference and finite element methods for numerical simulation of groundwater infiltration and pollution propagation. Šim˚unek et al. (2009) applied the finite element method to numerically solve the Richards equation and developed the commercial software Hydrus-1D. Zambra et al. (2011) constructed a finite volume method with high accuracy in space and time for solving nonlinear Richards equations. Ku et al. (2018) linearized the Richards' equation and developed a numerical solution to the unsteady groundwater infiltration issue by the collocation Trefftz method. Zeng proposed a modified Richards' equation and used FEM to more efficiently solve the variable saturated infiltration problem in heterogeneous soils. Chávez-Negrete et al. (2018) proposed an improved FDM combined with an adaptive stepsize Crank–Nicolson method for solving the Richards' equation. Li et al. (2022) established the convergence selection criteria of grid size and time step in finite element simulation of unsaturated infiltration, which can better enhance the numerical accuracy. Of course, other advanced numerical methods have certain advantages in terms of computational accuracy, computational efficiency or ease of implementation under certain conditions when solving variable saturated infiltration issues. These methods include hybrid finite element methods (Bergamaschi and Putti 1999), cellular automata methods (Mendicino et al. 2006), FDM with non-orthogonal grids (An et al. 2010), finite analysis methods (Zhang et al. 2016), meshless methods (Herrera et al. 2009; Ku et al. 2021), and Chebyshev spectral methods (Wu et al. 2020b). Table 3.1 lists the development of some classical numerical methods for solving Richards' equation.

In this chapter, numerical solutions of water infiltration equation in unsaturated soils are examined. Some improved methods such as fractional-order and Chebyshev spectral methods will be applied to investigate nonlinear infiltration in homogeneous and layered soils. The results will be compared with software simulation.

### **3.2 Numerical Solutions**

### *3.2.1 Finite Difference Method*

The FDM to solve Richards' equation first divides the solution area into differential grids, replaces the continuous solution domain with a finite number of grid nodes,


**Table 3.1** Some studies on the numerical solution of Richards' equation

(continued)


**Table 3.1** (continued)

and stores the variables to be solved (pressure head, water content, etc.) in each grid node. The differential term in the Richards' equation is replaced by the corresponding difference quotient, so that the partial differential equation is converted into an algebraic difference equation, and a differential equation system containing a finite number of unknown variables at discrete points is obtained, that is, a linear equation system. The solution of the linear equation system is obtained, and the numerical solution of the variables on the grid nodes is obtained. The principle is simple and easy to implement, and it is widely used. In Fig. 3.1, the distance along the *z*-axis is divided into *N* equal parts, the uniform grid step size is Δ*z*, and the spatial derivative first-order central difference scheme can be expressed as:

$$\left(\frac{\partial h}{\partial z}\right)\_i = \frac{h^\*\_{\ i+1} - h^\*\_{\ i-1}}{2\Delta z} \tag{3.1}$$

where *i* represents the discrete grid nodes along the *z*-axis.

The spatial derivative second-order central difference format is expressed as:

### 3.2 Numerical Solutions 45

$$\left(\frac{\partial^2 h}{\partial z^2}\right)\_i = \frac{h^\*|\_{i+1} - 2h^\*|\_i + h^\*|\_{i-1}}{\Delta z^2} \tag{3.2}$$

Additionally, the first-order backward difference format of the time derivative is given by:

$$\left(\frac{\partial h}{\partial t}\right)\_j = \frac{h^{\*j} - h^{\*j-1}}{\Delta t} \tag{3.3}$$

where Δ*t* is the discrete time step, and *j* is the time node.

Substituting Eqs. (3.1)–(3.3) into the Richards' equation, one can obtain (Zhu et al. 2022a):

$$\left(\frac{h\_{i+1}^{\*j} - 2h\_i^{\*j} + h\_{i-1}^{\*j}}{\Delta z^2}\right) + \alpha \left(\frac{h\_{i+1}^{\*j} - h\_{i-1}^{\*j}}{2\Delta z}\right) = c \left(\frac{h\_i^{\*j} - h\_i^{\*j-1}}{\Delta t}\right) \quad 1 \le i \le N - 1 \tag{3.4}$$

When the steady-state infiltration is incorporated, its finite difference format can be simplified to:

$$
\left(\frac{h\_{i+1}^\* - 2h\_i^\* + h\_{i-1}^\*}{\Delta z^2}\right) + \alpha \left(\frac{h\_{i+1}^\* - h\_{i-1}^\*}{2\Delta z}\right) = 0\tag{3.5}
$$

The system of linear equations formed by Eqs. (3.4)–(3.5) can be further written in matrix form:

$$\mathbf{A}\mathbf{h}^\* = \mathbf{b} \tag{3.6}$$

where **A** is a tridiagonal matrix of order (*N* − 1) × (*N* − 1); **h**<sup>∗</sup> and **b** are both column vectors of order (*N* − 1) × 1, and the first and last elements of vector **b** already contain boundary conditions. Equation (3.6) can be solved by relevant linear iterative methods, such as trigonometric decomposition method, Jacobi iterative method, Gauss–Seidel iterative method and relaxed iterative method.

### *3.2.2 Finite Volume Method*

The finite volume method (FVM), also known as the control volume method, the basic idea is to divide the calculation area into a series of non-repetitive control volumes, and make a control volume around each grid point, and the differential equation is divided into a control volume. For each control volume integral, a system of linear equations to be solved is obtained. For the one-dimensional uniform grid coordinates in Fig. 3.1, the control volume here can be assumed to be Δ*z*×1×1, that is, the distance in the other directions except the *z*-axis is unit length. Furthermore, by integrating the Richards' equation over the uniform control volume, one can obtain:

$$\int\_{i}^{i+\Delta z} \int\_{i}^{i+\Delta t} \frac{\partial \theta}{\partial t} \mathrm{d}t \,\mathrm{d}z = \int\_{i}^{i+\Delta t} \int\_{i}^{i+\Delta z} \frac{\partial}{\partial z} \left[ K\_{z}(h) \left( \frac{\partial h}{\partial z} + 1 \right) \right] \mathrm{d}z \,\mathrm{d}t \tag{3.7}$$

Furthermore, the discrete equation is written as:

$$\begin{split} \frac{K\_{i+1/2}^{j}(h\_{i+1}^{j} - h\_i^{j})}{\Delta z} - \frac{K\_{i-1/2}^{j}(h\_i^{j} - h\_{i-1}^{j})}{\Delta z} + K\_{i+1/2}^{j} - K\_{i-1/2}^{j} \\ = \Delta z C\_i^{j-1/2} \frac{h\_i^{j} - h\_i^{j-1}}{\Delta t} \end{split} \tag{3.8}$$

where *Ki*+1/2 and *Ki*−1/2 are the harmonic mean values of the permeability coefficients of adjacent nodes:

$$K\_{i+1/2} = \frac{2K\_iK\_{i+1}}{K\_i + K\_{i+1}} \tag{3.9}$$

$$K\_{i-1/2} = \frac{2K\_iK\_{i-1}}{K\_i + K\_{i-1}}\tag{3.10}$$

Equation (3.8) can be further simplified into the following matrix form:

$$\mathbf{A}(\mathbf{h}^j)\mathbf{h}^j = \mathbf{b}(\mathbf{h}^j, \mathbf{h}^{j-1})\tag{3.11}$$

For the solution of Eq. (3.8), the Picard iteration method is usually used to solve it.

### *3.2.3 Finite Element Method*

The finite element method (FEM) divides the solution domain into different elements, including triangular elements, rectangular elements, and quadrilateral elements. The equivalent integral equation of the problem can be obtained based on the variational principle and the orthogonalization principle of the weight function, and then the corresponding linear equation system can be obtained and solved iteratively. Usually, the Richards' equation solved by FEM is more than two-dimensional, then the twodimensional linearized Richards' equation considering the *x* and *z* directions can be written as (Zhu et al. 2022b, c):

$$\frac{\partial^2 h^\*}{\partial x^2} + \frac{\partial^2 h^\*}{\partial z^2} + \alpha \frac{\partial h^\*}{\partial z} = c \frac{\partial h^\*}{\partial t} \tag{3.12}$$

**Fig. 3.2** 2D infiltration model and rectangular element grid

As shown in Fig. 3.2, a rectangular element is used to divide the solution area. The number of discrete nodes along the *x* and *z* axes is *N*, then the discrete equation can be expressed in matrix form as follows:

$$\mathbf{G}\mathbf{h}^\* = -\mathbf{P}\frac{\mathbf{h}^{\*j} - \mathbf{h}^{\*j-1}}{\Delta t} \tag{3.13}$$

where **G** represents the total conduction matrix on the left side of Eq. (3.12); **P**  represents the storage matrix on the right side of Eq. (3.12); **h**<sup>∗</sup> is the array of the node pressure heads. Among them, the rectangular element (Fig. 3.2) composed of nodes *i*, *j*, *m*, and *n* only contributes to the *L*th row of **G** (*L* = *i*, *j*, *m*, and *n*). The *L*th row elements of matrix **G** in a rectangular element can be given by:

$$\mathcal{G}\_{L,i}^{\epsilon} = \int\_{-a}^{a} \int\_{-b}^{b} \left( \frac{\partial N\_i^{\epsilon}}{\partial x} \frac{\partial N\_L^{\epsilon}}{\partial x} + \frac{\partial N\_i^{\epsilon}}{\partial z} \frac{\partial N\_L^{\epsilon}}{\partial z} - \alpha\_{\text{g}} N\_i^{\epsilon} N\_L^{\epsilon} \right) \text{d}x \text{d}z \tag{3.14a}$$

$$G\_{L,j}^{\epsilon} = \int\_{-a}^{a} \int\_{-b}^{b} \left( \frac{\partial N\_j^{\epsilon}}{\partial x} \frac{\partial N\_L^{\epsilon}}{\partial x} + \frac{\partial N\_j^{\epsilon}}{\partial z} \frac{\partial N\_L^{\epsilon}}{\partial z} - \alpha\_g N\_j^{\epsilon} N\_L^{\epsilon} \right) \text{d}x \text{d}z \tag{3.14b}$$

$$G\_{L,m}^{\epsilon} = \int\_{-a}^{a} \int\_{-b}^{b} \left( \frac{\partial N\_m^{\epsilon}}{\partial x} \frac{\partial N\_L^{\epsilon}}{\partial x} + \frac{\partial N\_m^{\epsilon}}{\partial z} \frac{\partial N\_L^{\epsilon}}{\partial z} - \alpha\_s N\_m^{\epsilon} N\_L^{\epsilon} \right) \mathrm{d}x \mathrm{d}z \tag{3.14c}$$

48 3 Numerical Solutions to Infiltration Equation

$$\mathcal{G}\_{L,n}^{\epsilon} = \int\_{-a}^{a} \int\_{-b}^{b} \left( \frac{\partial N\_n^{\epsilon}}{\partial x} \frac{\partial N\_L^{\epsilon}}{\partial x} + \frac{\partial N\_n^{\epsilon}}{\partial z} \frac{\partial N\_L^{\epsilon}}{\partial z} - \alpha\_g N\_n^{\epsilon} N\_L^{\epsilon} \right) \mathrm{d}x \mathrm{d}z \tag{3.14d}$$

where *N <sup>e</sup> <sup>L</sup>*represents the element basis function; and *a* and *b* are half the length and width of the rectangular element, respectively. *N <sup>e</sup> <sup>L</sup>*can be formatted as:

$$N\_i^{\epsilon}(\mathbf{x}, z) = \frac{1}{4} \left( 1 - \frac{\mathbf{x}}{b} \right) \left( 1 - \frac{z}{a} \right) \tag{3.15a}$$

$$N\_j^{\epsilon}(\mathbf{x}, z) = \frac{1}{4} \left( 1 + \frac{\chi}{b} \right) \left( 1 - \frac{z}{a} \right) \tag{3.15b}$$

$$N\_m^{\varepsilon}(\mathbf{x}, z) = \frac{1}{4} \left( 1 + \frac{\mathbf{x}}{b} \right) \left( 1 + \frac{z}{a} \right) \tag{3.15c}$$

$$N\_n^{\epsilon}(\mathbf{x}, z) = \frac{1}{4} \left( 1 - \frac{\mathbf{x}}{b} \right) \left( 1 + \frac{z}{a} \right) \tag{3.15d}$$

Similarly, the *L*th row elements of **P** can also be formed as:

$$P\_{L,i}^{\varepsilon} = c \int\_{-a}^{a} \int\_{-b}^{b} N\_i^{\varepsilon} N\_L^{\varepsilon} \mathrm{d}x \mathrm{d}z \tag{3.16a}$$

$$P\_{L,j}^{\epsilon} = c \int\_{-a}^{a} \int\_{-b}^{b} N\_j^{\epsilon} N\_L^{\epsilon} \mathrm{d}x \mathrm{d}z \tag{3.16b}$$

$$P\_{L,m}^{\epsilon} = c \int\_{-a}^{a} \int\_{-b}^{b} N\_m^{\epsilon} N\_L^{\epsilon} \mathrm{d}x \mathrm{d}z \tag{3.16c}$$

$$P\_{L,n}^{\varepsilon} = c \int\_{-a}^{a} \int\_{-b}^{b} N\_n^{\varepsilon} N\_L^{\varepsilon} \mathrm{d}x \mathrm{d}z \tag{3.16d}$$

Furthermore, the contributions of all rectangular elements are summed and assembled to form matrices **G** and **P**. Generally, Eqs. (3.14)–(3.16) are solved using Gaussian integration method (Wang and Anderson 1982). The Gaussian integration method makes the definite integral equal to a weighted sum over a finite number of points, and for a quadratic polynomial, the general expression for the double integral is as follows:

3.2 Numerical Solutions 49

$$\int\limits\_{-1}^{1}\int\limits\_{-1}^{1}f(\varepsilon,\eta)\mathrm{d}\varepsilon\mathrm{d}\eta = f(\varepsilon\_{1},\eta\_{1}) + f(\varepsilon\_{2},\eta\_{1}) + f(\varepsilon\_{1},\eta\_{2}) + f(\varepsilon\_{2},\eta\_{2})\tag{3.17}$$

Among them, the four Gauss points are determined by ε1 = − 1/ <sup>√</sup>3, ε2 <sup>=</sup> <sup>1</sup>/ <sup>√</sup>3, η1 = − 1/ <sup>√</sup>3, η2 <sup>=</sup> <sup>1</sup>/ <sup>√</sup>3, and the ownership coefficients are all equal to 1. The integral of Eqs. (3.14)–(3.16) can be transformed into the form of Eq. (3.17) by the following variable transformation:

$$
\varepsilon = \frac{x}{b} \text{ and } \eta = \frac{z}{a} \tag{3.18}
$$

Then, d*z* = *a*dε, d*x* = *b*dη, the integral limit is from−1 to 1, and the interpolation function (Eq. 3.15) in the coordinates (ε, η) is transformed into:

$$\overline{N}\_i^{\varepsilon}(\varepsilon,\eta) = \frac{1}{4}(1-\varepsilon)(1-\eta) \tag{3.19a}$$

$$\overline{N}\_j^{\varepsilon}(\varepsilon,\eta) = \frac{1}{4}(1+\varepsilon)(1-\eta) \tag{3.19b}$$

$$\overline{N}\_{\mathfrak{m}}^{\varepsilon}(\varepsilon,\eta) = \frac{1}{4}(1+\varepsilon)(1+\eta) \tag{3.19c}$$

$$\overline{N}\_n^{\epsilon}(\varepsilon,\eta) = \frac{1}{4}(1-\varepsilon)(1+\eta) \tag{3.19d}$$

In summary, the integral transformation of Eqs. (3.14)–(3.16) is described as:

$$\int\_{-a}^{a} \int\_{-b}^{b} f(\mathbf{x}, \mathbf{y}) \mathrm{d}\mathbf{x} \mathrm{d}\mathbf{y} = ab \int\_{-1}^{1} \int\_{-1}^{1} f(\varepsilon, \eta) \mathrm{d}\varepsilon \mathrm{d}\eta \tag{3.20}$$

where *f* (ε, η) is the transformation of *f* (*x*,*z*) at the coordinate (ε, η).

### *3.2.4 Numerical Approximation to the Fractional-Time Richards' Equation*

Fractional Richards' equations can be further divided into space fractional, time fractional, and space–time fractional partial differential equations. The time fractional Richards' equation is expressed as follows:

$$C(h)\frac{\partial^{\mathcal{V}}h}{\partial t^{\mathcal{Y}}} = \frac{\partial}{\partial z}\Big[K(h)\Big(\frac{\partial h}{\partial z} + 1\Big)\Big] \tag{3.21}$$

where γ is the order of the time derivative. When fractional order γ = 1, Eq. (3.21) degenerates into the classical Richard equation.

The Caputo fractional derivative with respect to the function *f* (*x*) can be defined as (Pachepsky et al. 2003):

$$\, \_0^C D\_x^{\gamma} f(\mathbf{x}) = \frac{1}{\Gamma(n - \gamma)} \int\_0^\mathbf{x} \frac{f^{(n)}(\mathbf{r})}{(\mathbf{x} - \mathbf{r})^{\gamma - n + 1}} d\mathbf{r} \tag{3.22}$$

where *<sup>C</sup>* 0 *D*<sup>γ</sup> *<sup>x</sup>*represents the γ -order Caputo fractional derivative; ┌ is the gamma function, and *n* is a positive integer (representing an integer derivative). In addition, the Riemann–Liouville (R-L) fractional derivative can be defined as (Su 2014):

$$\, \_0^{RL} D\_x^{\gamma} f(\mathbf{x}) = \frac{1}{\Gamma(n - \gamma)} \frac{\mathbf{d}^n}{\mathbf{d}x^n} \int\_0^\mathbf{x} \frac{f(\mathbf{r})}{(\mathbf{x} - \mathbf{r})^{\gamma - n + 1}} \mathbf{d}\mathbf{r} \tag{3.23}$$

where *RL*  0 *D*<sup>γ</sup> *<sup>x</sup>*represents the γ -order R-L fractional derivative.

It can be seen from Eqs. (3.22)–(3.23) that in the definition of R-L fractional derivative, the fractional integral is first obtained and then the integer derivative is obtained, while the definition of Caputo derivative is to first obtain the integer derivative, and then the fractional integral is calculated. There is a great connection between the two, and there are also essential differences, which are mainly reflected in the actual physical models. In the process of actually solving the initial value problem of differential equations, Caputo derivatives are more widely used and have more physical background than R-L derivatives.

For the solution of Eq. (3.21), the finite difference method is used for numerical discretization, and meshing is conducted. Let η and Δ*z* be the time and space steps, respectively. The simulation time is divided into *M* equal parts, and the *z*-axis is divided into *N* equal parts:

$$t\_m = m\eta, \quad m = 0, 1, \dots, M \tag{3.24}$$

$$z\_i = i \,\Delta z, \quad i = 0, 1, \ldots N \tag{3.25}$$

For the left-hand fractional derivative term of Eq. (3.21), the Caputo fractional derivative is defined as:

$$\frac{\partial^{\mathcal{V}}h(z,t)}{\partial t^{\mathcal{V}}} = \frac{1}{\Gamma(1-\mathcal{V})} \int\_0^t \frac{\partial h(z,\tau)}{\partial \tau} \frac{d\tau}{(t-\tau)^{\mathcal{V}}}\tag{3.26}$$

The integral of the pressure head derivative in Eq. (3.26) can be directly approximated by the numerical differential equation, which is deduced as follows:

$$\frac{\partial^{\mathcal{V}}h(z,t\_m)}{\partial t^{\mathcal{V}}} = \frac{1}{\Gamma(1-\mathcal{V})} \int\_0^{t\_m} \frac{h'(\tau)d\tau}{(t\_m-\tau)^{\mathcal{V}}} = \frac{1}{\Gamma(1-\mathcal{V})} \sum\_{j=0}^{m-1} \int\_{t\_j}^{t\_{j+1}} \frac{h'(t\_m-\tau)d\tau}{\tau^{\mathcal{V}}}$$

$$\approx \frac{1}{\Gamma(1-\mathcal{V})} \sum\_{j=0}^{m-1} \frac{h(t\_m-t\_j)-h(t\_m-t\_{j+1})}{\eta} \int\_{t\_j}^{t\_{j+1}} \tau^{-\mathcal{V}}d\tau$$

$$= \frac{\eta^{-\mathcal{V}}}{\Gamma(2-\mathcal{V})} \sum\_{j=0}^{m-1} \left(h\_{m-j}-h\_{m-j-1}\right) \left[(j+1)^{1-\mathcal{V}}-j^{1-\mathcal{V}}\right] \tag{3.27}$$

Equation (3.27) is further simplified as:

$$\frac{\partial^{\gamma}h(z,t\_m)}{\partial t^{\gamma}} = \frac{\eta^{-\gamma}}{\Gamma(2-\gamma)}\sum\_{j=0}^{m-1}b\_j^{(\gamma)}\left(h\_{m-j} - h\_{m-j-1}\right) \tag{3.28}$$

where *b*(γ ) *<sup>j</sup>*= (*j* + 1) <sup>1</sup>−<sup>γ</sup> <sup>−</sup> *<sup>j</sup>* <sup>1</sup>−<sup>γ</sup> .

For finite difference discretization on the right side of Eq. (3.21), one obtains:

$$\frac{\partial}{\partial z} \left[ K(h) \left( \frac{\partial h}{\partial z} + 1 \right) \right] = \frac{K\_{i+1/2}^{m+1} (h\_{i+1}^{m+1} - h\_i^{m+1})}{\Delta z^2} - \frac{K\_{i-1/2}^{m+1} (h\_i^{m+1} - h\_{i-1}^{m+1})}{\Delta z^2} $$

$$+ \frac{K\_{i+1/2}^{m+1} - K\_{i-1/2}^{m+1}}{\Delta z} \tag{3.29}$$

*Ki*+1/2 and *Ki*−1/2 can be expressed as:

$$K\_{i+1/2} = \frac{2K\_iK\_{i+1}}{K\_i + K\_{i+1}}\tag{3.30}$$

$$K\_{i-1/2} = \frac{2K\_iK\_{i-1}}{K\_i + K\_{i-1}}\tag{3.31}$$

Combining Eqs. (3.28) and (3.29), the discrete format of the fractional-time Richards' equation can be obtained as follows:

$$\begin{split} & \frac{\eta^{-\mathcal{V}}}{\Gamma(2-\mathcal{V})} \mathcal{C}(h)\_{i}^{m-1/2} \sum\_{j=0}^{m-1} b\_{j}^{(\mathcal{V})} \{h\_{m-j,i} - h\_{m-j-1,i}\} \\ &= \frac{K\_{i+1/2}^{m} (h\_{i+1}^{m} - h\_{i}^{m})}{\Delta z^{2}} - \frac{K\_{i-1/2}^{m} (h\_{i}^{m} - h\_{i-1}^{m})}{\Delta z^{2}} + \frac{K\_{i+1/2}^{m} - K\_{i-1/2}^{m}}{\Delta z} \end{split} \tag{3.32}$$

where *C*(*h*) *m*−1/2 *<sup>i</sup>* represents the average value of the specific moisture capacity (*C*(*h*) *m <sup>i</sup>*) at the previous time step and the specific moisture capacity (*C*(*h*) *m*−1,*k <sup>i</sup>* ) of the current iteration step (*k* ≥ 1) of the current time step.

Equation (3.32) can be simplified into matrix form **Ax** = **b**, and then iteratively solved by Picard method. In order to evaluate the fitting effect of the proposed time fractional model, two indicators are selected, namely the root mean square error (RMSE) and the relative error (RE):

$$\text{RSE} = \left[ \frac{1}{N - 1} \sum\_{i = 1}^{N - 1} \left( \frac{h\_i - h\_i^\*}{h\_i^\*} \right)^2 \right]^{0.5} \tag{3.33}$$

$$\text{RE} = 100 \times \left[ \frac{1}{N - 1} \sum\_{i = 1}^{N - 1} \left( \left| \frac{h\_i - h\_i^\*}{h\_i^\*} \right| \right) \right] \tag{3.34}$$

where *hi* is the numerical solution of Richards' equation, and *h*<sup>∗</sup> *i* is the exact solution of Richards' equation. The smaller the values of the two errors, the higher the computational accuracy of the proposed approach.

### **3.3 An Improved Method for Non-uniform Spatial Grid**

In the numerical solution of Richards' equation, FDM can be used for numerical discretization and iterative solution. However, to obtain a more reliable numerical solution, the space step size of conventional uniform grid (Fig. 3.1) is often small, particularly under some unfavorable numerical conditions, such as infiltration into dry and layered soils with greatly different permeability coefficients, this makes the calculation time-consuming and even the accuracy cannot be improved very well. The realization of unstructured space grids and dynamic grid methods in some studies is often complicated and inappropriate (Chávez-Negrete et al. 2018; Dolejší et al. 2019). Therefore, this chapter proposes an improved FDM numerical discretization process using a non-uniform Chebyshev space grid, which is compared with the traditional uniform space grid. This method can provide a certain reference for the numerical simulation of unsaturated infiltration.

(1) Uniform grid method

The Richards' equation is directly numerically discretized by the finite difference method of uniform grid, and one can gain (Zhu et al. 2020):

$$\frac{\left[K\_{i+1/2}^{j}(h\_{i+1}^{j}-h\_{i}^{j})-K\_{i-1/2}^{j}(h\_{i}^{j}-h\_{i-1}^{j})\right]}{\Delta z^{2}}+\frac{K\_{i+1/2}^{j}-K\_{i-1/2}^{j}}{\Delta z}$$

$$=\frac{C\_{i}^{j-1/2}\left(h\_{i}^{j}-h\_{i}^{j-1}\right)}{\Delta t}\tag{3.35}$$

Comparing Eq. (3.29), it can be found that the mathematical meaning of Eqs. (3.35) and (3.29) is consistent. When its steady-state infiltration is considered, the FDM format can be simplified as (Zhu et al. 2020):

$$\frac{\left[K\_{i+1/2}(h\_{i+1} - h\_i) - K\_{i-1/2}(h\_i - h\_{i-1})\right]}{\Delta z^2} + \frac{K\_{i+1/2} - K\_{i-1/2}}{\Delta z} = 0 \qquad (3.36)$$

### (2) Improved Chebyshev space grid method

Because the pressure head in the unsaturated infiltration problem often has a large change at the boundary and the soil layer interface, a finer mesh is usually required for densification. However, the uniform grid generates too many computing grid nodes in the process of densification, and the computation is time-consuming and even not accurate enough. Then, a Chebyshev grid coordinate is proposed as (Wu et al. 2020b; Zhu et al. 2020):

$$z\_i = \cos(i\pi/N) \times \frac{L}{2} + \frac{L}{2}, \quad i = N, N-1, \dots, 0\tag{3.37}$$

where *L* is the thickness of the soil layer.

In Fig. 3.1, the Chebyshev grid nodes are only highly refined at the interface, greatly reducing the number of grid nodes. Furthermore, the FDM discrete format combined with the Chebyshev grid can be expressed as:

$$\frac{1}{\delta z\_i} \left[ \frac{K\_{i+1/2}^j (h\_{i+1}^j - h\_i^j)}{\Delta z\_{i+1}} - \frac{K\_{i-1/2}^j (h\_i^j - h\_{i-1}^j)}{\Delta z\_i} \right] + \frac{K\_{i+1/2}^j - K\_{i-1/2}^j}{\delta z\_i}$$

$$= C\_i^{j-1/2} \frac{h\_i^j - h\_i^{j-1}}{\Delta t} \tag{3.38}$$

where Δ*zi* represents the unequal spacing between the Chebyshev grid nodes; δ*zi* is the unequal spacing of the computing nodes (Fig. 3.1).

Equation (3.38) can also be simplified into a matrix form such as Eq. (3.13). Due to the nonlinear relationship between permeability coefficient and water content, the coefficient matrix **A** needs to be repeatedly evaluated by nonlinear iteration method after numerical discretization. Among them, Picard method is a more classical and practical nonlinear iterative method. Based on the Chebyshev grid discretization format of the Richards' equation, a program for unsaturated infiltration was developed using the MATLAB (R2014a) language.

### *3.3.1 Validation Example*

This test describes one-dimensional transient unsaturated infiltration in homogeneous unsaturated soil (Ku et al. 2018; Zhu et al. 2022a), the soil thickness *L* = 10 m, exponential model parameters are: α = 1 × 10−4, θs = 0.50, θr = 0.11, and the saturated permeability coefficient *k*s = 2.5 × 10−8 m/s. Furthermore, the boundary conditions can be given by:

$$h(z=0) = h\_{\mathbb{d}}\tag{3.39}$$

$$h(z=10) = 0\tag{3.40}$$

To verify the computational accuracy of the proposed method, the following error indicators are used in the comparative analysis, namely the root mean square error (RSE), the relative error (RE), and maximum relative error (MRE):

$$\text{RSE} = \left[ \frac{1}{N - 1} \sum\_{i = 1}^{N - 1} \left( \frac{h\_i - h\_i^\*}{h\_i^\*} \right)^2 \right]^{0.5} \tag{3.41}$$

$$\text{RE} = 100 \times \left[ \frac{1}{N - 1} \sum\_{i = 1}^{N - 1} \left( \left| \frac{h\_i - h\_i^\*}{h\_i^\*} \right| \right) \right] \tag{3.42}$$

$$\text{MRE} = 100 \times \left[ \max \left| \frac{h\_i - h\_i^\*}{h\_i^\*} \right|\_{i=1}^{N-1} \right] \tag{3.43}$$

where *hi* is the numerical solution of Richards' equation and *h*<sup>∗</sup> *<sup>i</sup>*is the exact solution of Richards' equation. The smaller the values of the three errors, the higher the computational accuracy of the proposed method. The total simulation time was set to 5 h, the number of discrete nodes was set to 100, 150, and 200, and the time steps were set to 0.01, 0.02, and 0.04 h, respectively.

In Fig. 3.3, the numerical solutions are derived using two grid methods under the conditions of Δ*t* = 0.01 h and *N* = 200 and compared with the analytical solutions. The numerical solution obtained by the uniform grid method has a large deviation from the exact solution, particularly after *t* > 2 h (Fig. 3.3a), while the numerical solution obtained by the Chebyshev grid method is in good agreement with the analytical solution (Fig. 3.3b).

Figure 3.4a represents the maximum relative error (MRE) obtained by different grid methods at different time steps when the number of nodes *N* is 100. The MRE obtained by the Chebyshev grid method ranges from 0.6 to 3.5%, and it first decreases and then increases with the increase of time. When *t* < 4 h, the MRE decreases with decreasing time step. The MRE obtained by the uniform grid method ranges from 1.3 to 49%, particularly when *t* > 1 h, the MRE increases over time and is much larger than that obtained by the Chebyshev grid.

Figure 3.4b shows the MRE obtained by different grid methods under different numbers of discrete grid nodes when the time step Δ*t* = 0.01 h. It can be found that the MRE obtained by the uniform grid as time increases. In addition, the MRE of the two methods has a decreasing trend with the increase of the number of grid nodes *N*.

**Fig. 3.3** Comparison of numerical and analytical solutions obtained by different grid methods

In Table 3.2, the RSE of both methods and the RE of the uniform grid decrease with the increase of *N*, while the RE of the Chebyshev grid increases with increasing *N*. However, it can be seen from the numerical value that the RSE of the Chebyshev grid is nearly 100 times different from the uniform grid method, and the RE is more than 10 times different from the uniform grid method. This test indicates that the proposed Chebyshev grid method is not limited by the number of grid nodes

**Fig. 3.4** Comparison of the maximum relative error of different grid methods under different numerical conditions

to improve the accuracy, that is, the proposed method achieves higher numerical accuracy with fewer discrete nodes, and has a smaller computational cost.

### *3.3.2 Unsaturated Infiltration in Layered Soils*

The mathematical model is shown in Fig. 3.5. The model parameters of the two-layer soils are set to θs = 0.35, θr = 0.14, α = 8 × 10−3. The thickness of soil layer 1 and


**Table 3.2** Numerical accuracy at *t* = 5h

soil layer 2 are both set to 5 m, the number of discrete nodes in each layer is 40, and the boundary conditions are consistent with Sect. 3.3.1, where *h*d = − 103 m.

In layered soils, due to the influence of different soil types, different combinations have a great influence on unsaturated infiltration. Table 3.3 lists the saturated permeability coefficients in different saturated soils. In order to further verify the applicability of the proposed grid method, it is assumed that the saturated permeability coefficient of soil layer 1 is 10−1 m/s, and the saturated permeability coefficient of the second layer gradually changes from coarse sandy soil to clay in Table 3.3, the saturated permeability coefficient of the second layer ranges from 10−2 to 10−9 m/s.


It can be seen from Fig. 3.6a that the numerical solutions obtained by the uniform grid and the Chebyshev grid are approximate, and both can well simulate the unsaturated infiltration process in case 2. Figure 3.6b depicts the numerical results of case 4. The ratio of the saturated permeability coefficients of the upper and lower soil layers is in the order of 104. The uniform grid cannot accurately represent the pressure head at the interface, but the numerical solution obtained by the Chebyshev grid can accurately describe the pressure head variation at the interface. In Fig. 3.6c, the Chebyshev grid can also accurately describe the variations in the pressure head at the interface under case 6 (Table 3.4).

Similar to the two-layer soil, the proposed improved grid method is applied to the three-layer unsaturated soils. The mathematical model is described in Fig. 3.7, where *L*1 = *L*3 = 4 m and *L*2 = 2 m. The number of discrete nodes in each layer of soil is 40, the parameter θs is set to 0.46, the three-layer unsaturated soil permeability coefficient is listed in Table 3.5, and the boundary conditions can be expressed as: *h*(*z* = 0) = 0, *h*(*z* = 10 m) = − 1000 m. Other numerical conditions are consistent with the two-layered soils.

In Fig. 3.8, the proposed Chebyshev grid method can better characterize the pressure head change between the two interfaces in case 11 compared to the uniform grid method. This test further demonstrates that the proposed Chebyshev grid method can obtain more reliable numerical solutions with fewer grid nodes under some unfavorable infiltration conditions, particularly when the saturated permeability coefficient changes greatly in layered soils.

### **3.4 Application of Chebyshev Spectral Method to Richards' Equation**

### *3.4.1 Chebyshev Spectral Method*

The Chebyshev spectral method (CSM) was developed by Gottlieb et al. (1978). It is widely used to solve numerical problems expressed by partial differential equations. The advantage of the CSM is the use of non-uniform mesh discretization and the Chebyshev differentiation matrix to improve the accuracy of the numerical simulation. In Fig. 1.2, the *N* + 1 Chebyshev point coordinates in the interval [− 1, 1] can be expressed as:

$$z\_j = \cos(j\pi/N), \quad j = 0, 1, \ldots N \tag{3.44}$$

The first-order derivative of interpolation function *p* (**z**) can be written as:

$$p'(\mathbf{z}) = \mathbf{D}\_N \mathbf{h} \tag{3.45}$$


**Table 3.4** Permeability coefficient for cases 1–8

**Fig. 3.7** One-dimensional infiltration model of three-layer unsaturated soil


**Table 3.5** Permeability coefficient for cases 9–13

where **z** is represented as the *N* + 1-dimensional vector (*z*0,*z*1,...,*zN* ) T , **h** represents the *N* + 1 dimensional vector (*h*0, *h*1,..., *hN* ) T composed of the pressure head, and **D***N* is an (*N* + 1) × (*N* + 1) Chebyshev differentiation matrix.

The expression of each element in Chebyshev differentiation matrix **D***N* for any number of nodes (*N* + 1) is given by:

$$(\mathbf{D}\_N)\_{00} = \frac{2N^2 + 1}{6}, \quad (\mathbf{D}\_N)\_{NN} = -\frac{2N^2 + 1}{6}, \tag{3.46a}$$

**Fig. 3.8** Comparison of numerical solutions for case 11

$$(\mathbf{D}\_N)\_{jj} = \frac{-z\_j}{2\left(1 - z\_j^2\right)}, \quad j = 1, \ldots, N - 1 \tag{3.46b}$$

$$(\mathbf{D}\_N)\_{ij} = \frac{c\_i}{c\_j} \frac{(-1)^{i+j}}{z\_i - z\_j}, \quad (i \neq j)i, j = 0, \dots, N,\tag{3.46c}$$

where (**D***N* )*ij* represents the element in row *i* + 1 and column *j* + 1 of Chebyshev differentiation matrix **D***N* for the first derivative and,

$$c\_i = \begin{cases} 2, \ i = 0, N \\ 1, \ i = 1, \ldots, N - 1 \end{cases} \tag{3.47}$$

To obtain the Chebyshev differentiation matrix for the second derivative, the matrix can be directly squared to **D***N* . Similarly, Chebyshev differentiation matrices for higher-order derivatives can be obtained as follows:

$$\frac{\partial}{\partial z} \to \mathbf{D}\_N \tag{3.48a}$$

$$\frac{\partial^2}{\partial z^2} \to \mathbf{D}\_N^2 \tag{3.48b}$$

$$\frac{\partial^3}{\partial z^3} \rightarrow \mathbf{D}\_N^3 \tag{3.48c}$$

Because **D***N* is obtained in the interval [− 1, 1], it is necessary to scale **D***N* to other intervals, that is,

$$\frac{\partial}{\partial z} \rightarrow \frac{2}{L} \mathbf{D}\_N \tag{3.49a}$$

$$\frac{\partial^n}{\partial z^n} \to \left(\frac{2}{L}\mathbf{D}\_N\right)^n\tag{3.49b}$$

Through the above process, the partial differential equation including the RE can be easily and efficiently solved. Thus, the governing equation (Eq. 2.1) can be rewritten as:

$$\frac{1}{c}\mathbf{D}\_N^2 \mathbf{h}^\* + \frac{\alpha\_g}{c} \mathbf{D}\_N \mathbf{h}^\* = \frac{\partial \mathbf{h}^\*}{\partial t} \tag{3.50}$$

In Eq. (3.50), Eq. (2.1) has been transformed into an ordinary differential equation (ODE) by introducing Chebyshev differentiation matrices. In this book, the variablestep four- and fifth-order Runge–Kutta method is used to solve Eq. (3.50), which is the ODE solver in MATLAB. Meanwhile, the 2D linear RE can also be solved using the same method.

To evaluate the performance of the proposed method, an *L*2 norm is used to indicate the error of solution *h*∗:

$$L\_2(h^\*) = \left\| h^{\*,a}(z,t) - h^\*(z,t) \right\|\_2 \tag{3.51}$$

where *h*∗,*<sup>a</sup>*(*z*, *t*) and *h*∗(*z*, *t*) are one-dimensional (1D) analytical and numerical solutions, respectively. Simultaneously, the numerical and analytical solutions with spatial and temporal grids can produce an absolute error defined as follows:

$$\text{err}\_A(h^\*) = h^{\*,a}(z,t) - h^\*(z,t). \tag{3.52}$$

### *3.4.2 Validation Example*

Test 1 represents transient infiltration in homogeneous unsaturated soils. In the mathematical model, the soil thickness (*L*) is assumed to be 10 m. The model parameters are θs = 0.50, θr = 0.11, α = 1 × 10−4, and *k*s = 9 × 10−5 m/h (Zhu et al. 2019). The governing equation is described as Eq. (2.1).

When water infiltrates into a soil mass, ponding on the ground maintains the pressure head at zero (Green and Ampt 1911). The upper and lower boundary conditions can be expressed as:

$$h(z=0,t) = h\_{\mathrm{d}}\tag{3.53}$$

$$h(z=L,t) = 0\tag{3.54}$$

The effect of Δ*z* and *N* on the results of Test 1 was analyzed. The time step was set to 0.01 h. Different numbers of nodes or grid sizes, that is, *N*/Δ*z* = 20/0.5 (m), *N*/Δ*z* = 40/0.25 (m), and *N*/Δ*z* = 80/0.125 (m) were selected to investigate the influence of the different methods on *L*2(*h*∗). The total simulation time was 5 h. Simultaneously, three different initial conditions *h*0 were chosen to evaluate their influence on the numerical accuracy and stability of the FDM and CSM.

In Fig. 3.9, *L*2(*h*∗) of the CSM were consistently smaller than that of the conventional FDM. The numerical accuracy increased as the mesh was refined, and the accuracy of the traditional FDM was significantly affected by the initial conditions. However, the proposed method (CSM) was less affected, and the numerical accuracy was stable in the order of 10−6–10−7 in this example. The result indicates that the CSM was more robustness than the FDM. Figure 3.10 demonstrates the comparison of the computed pressure head profiles with different grid sizes at *t* = 2 h. Compared with the analytical solution, it can be easily seen that the numerical accuracy of the CSM was less influenced by the grid sizes than that of the FDM (Fig. 3.10). That is, the CSM achieved better numerical results with fewer mesh nodes.

The pressure head profiles were obtained from the computed results over time under *h*0 = − 100,000 m (Fig. 3.11). The number of nodes or grid sizes was 40/0.25 (m). Figure 3.11 illustrates that the numerical solutions of the CSM agree well with the analytical solutions over time, whereas the numerical solutions of the FDM have larger errors than the analytical solutions. In Fig. 3.12, the minimum absolute errors of the proposed CSM and the FDM are approximately 10−11 and 10−6, respectively. Consequently, the extended CSM is characterized by higher accuracy than traditional FDM with fewer nodes.

### **3.5 Conclusions**

This chapter first summarizes the commonly used spatial numerical discrete methods of Richards' equation, including FDM, FVM, and FEM. Furthermore, this chapter proposes an improved FDM numerical discretization process using a non-uniform Chebyshev space grid and compares it with the numerical results and analytical solutions obtained from a conventional uniform space grid. Additionally, based on the non-uniform Chebyshev grid method, the Chebyshev spectral method is proposed to solve the Richards' equation, and the following conclusions can be obtained:

(1) The proposed Chebyshev grid method, while keeping the number of discrete nodes unchanged, uses a cosine function to ingeniously densify the interfaces on both sides and then numerically solve it to obtain a more reliable numerical solution. The numerical results indicate that the numerical solutions obtained

**Fig. 3.9** Comparison of the calculation accuracy for solving 1D transient infiltration with different grid sizes and initial conditions using different methods: **a** *h*<sup>0</sup> = − 10 m; **b** *h*0 = − 1000 m; **c** *h*0 = − 100,000 m

**Fig. 3.11** Comparison of the calculation results (pressure head) from different methods when *h*<sup>0</sup> = − 100,000 m

**Fig. 3.12** Absolute error in the results computed using the CSM and FDM relative to the analytical solution for the 1D transient infiltration at different simulated times

by this method are in good agreement with the analytical solutions under some unfavorable numerical conditions, such as infiltration into a dry soil. At the same time, the results illustrate that the proposed Chebyshev grid method can obtain higher numerical accuracy with a smaller number of nodes than the conventional uniform grid, and the computational cost is small.

(2) The Chebyshev spectral method (CSM) based on the discretization of Chebyshev grid is extended to simulate unsaturated infiltration in porous media. Chebyshev differential matrix can construct *n*-order Chebyshev differential matrix conveniently and quickly. This is the first time that the CSM has been used successfully to solve 1D and 2D transient infiltration problems in unsaturated soils. Compared with the traditional FDM, the CSM was more efficient, and it achieved higher accuracy with a lower-resolution grid. The CSM was not sensitive to the initial conditions.

### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Chapter 4 Improved Linear and Nonlinear Iterative Methods for Rainfall Infiltration Simulation**

71

### **4.1 Introduction**

The linear infiltration equations obtained by discretizing Richards' equation need to be solved iteratively, including two approaches of linear and nonlinear iterations. The first method is to use numerical methods to directly numerically discretize Richards' equations to obtain nonlinear ordinary differential equations and then use nonlinear iterative methods to iteratively solve, such as Newton's method (Radu et al. 2006), Picard method (Lehmann and Ackerer 1998), and the *L*-method (List and Radu 2016). Among them, Paniconi and Putti (1994) found that Newton's method did not converge systematically, and sometimes the numerical results could be completely wrong. Casulli and Zanolli (2010) proposed a nested Newton method that can obtain quadratic convergence rates for arbitrary discrete time steps and all flow regimes. Brenner and Cancès (2017) introduced a new parameter to modify the Richards equation, which will be more robust when solving the equation using Newton's method. The Picard method can be considered as a simplified Newton method, which linearly converges. Since Newton's method increases algebraic complexity and computational cost for assembling the derivative terms in the Jacobian matrix (Zeng et al. 2018), the modified Picard method (Celia et al. 1990) has been widely studied for *h*-based and mixed forms of Richards' equation. However, since the coefficient matrix in each iteration needs to be re-evaluate in the Picard method, there is computationally expensive (Farthing and Ogden 2017). Some unfavorable numerical conditions such as water infiltration into dry and layered soils with very different hydraulic conductivities lead to numerical difficulties or even non-convergence in the conventional Picard method (Zha et al. 2017). Therefore, there has been increasing attention on the improvement of numerical discretization and computational efficiency of the Picard method (Celia et al. 1990; Zhu et al. 2020).

The second approach first needs to convert the Richards' equation into a linearized partial differential equation (Wu et al. 2020a, b), and then numerical methods are utilized to perform numerical discretization to obtain a set of linear equations and

<sup>©</sup> The Author(s) 2023

L. Wu and J. Zhou, *Rainfall Infiltration in Unsaturated Soil Slope Failure*, SpringerBriefs in Applied Sciences and Technology, https://doi.org/10.1007/978-981-19-9737-2\_4

solve them iteratively. The classical linear iterative methods can efficiently solve systems of linear equations, such as the classic Jacobi iterative method (Hagemam and Young 1981), Gauss–Seidel iterative method (GS), and successive overrelaxation method (SOR). The GS method is an improvement on the Jacobi method, which usually takes less computation time, but it still requires more iterations to obtain a numerical solution than the SOR method. The Chebyshev semi-iterative method (CSIM) can improve the iterative convergence rate by introducing Chebyshev polynomials (Arioli and Scott 2014), but the computational cost of CSIM in each iteration is lower than that of traditional GS. Additionally, the Krylov subspace iterative method is also a common iterative method for solving linear equations, such as the conjugate gradient method, the generalized minimum residual method (GMRES), and the modified bi-conjugate gradient method (Dehghan and Mohammadi-Arani 2016). Theoretical analysis indicated that Hermitian/skew-Hermitian splitting iterative approaches (HSS) absolutely converge to a unique solution to a system of linear equations (Bai et al. 2010). HSS is an efficient linear iterative method for solving sparse non-Hermitian positive-definite equations (Dehghan and Shirilord 2019). Although HSS can converge unconditionally, it is time-consuming and impractical in real-world computation. Therefore, Bai et al. (2004) proposed the inexact Hermitian/skew-Hermitian splitting iterative approaches (IHSS) combined with the Krylov subspace iterative method to improve computational efficiency. However, IHSS is still time-consuming to solve linear equations for finer mesh discretization. Due to some shortcomings of conventional iterative methods, it is of great significance to improve numerical stability, accuracy, and convergence rate (Lott et al. 2012; Deng and Wang 2017; Mitra and Pop 2019; Illiano et al. 2021; Su et al. 2022).

The improvement methods for common accelerated convergence include extrapolation, error correction, Chebyshev polynomial acceleration (Arioli and Scott 2014), and preconditioning technique (Benzi 2002). For example, Lott et al. (2012) studied the effectiveness of the improved Picard method by the Anderson acceleration method for the variable saturated flow problem. Arioli and Scott (2014) analyzed an improved iterative method of Chebyshev polynomials that can be applied to accelerate iterative optimization without losing numerical stability. Wang and Zhang (2003) proposed a class of parallel multi-step sequential preprocessing strategies, which can significantly improve the computational efficiency and stability of solving linear equations. Briggs et al. (2000) proposed a multi-grid correction method to quickly eliminate the iterative error in the iterative process, thereby obtaining a faster convergence rate. Recently, preconditioning methods can effectively reduce the condition number of the iterative matrix of linear equations, thereby improving the computational convergence rate (Benzi 2002; Zhu et al. 2022a, b, c), including the left preconditioning, right preconditioning, and two-side preconditioning (Liu et al. 2015). The preconditioning technique converts the original system of linear equations into a system that is easier to solve, thereby improving the convergence rate of the iterative process (Benzi 2002).

With the higher accuracy of the numerical solutions to solve the infiltration problem, the scale of the system of algebraic equations increases, and then the calculation time also increases. Therefore, the numerical methods and iterative methods are constantly improving and innovating. The research based on traditional methods and newly developed and new algorithms is of great significance for the crossdevelopment of geotechnical engineering, geological engineering, and other majors. In this chapter, the Richards' equation can be transformed into a linearized Richards' equation using an exponential function. Furthermore, the numerical method is used to discretize the linear RE, and the linear and nonlinear iterative methods are employed to evaluate computational efficiency and accuracy using different iterative methods.

### **4.2 A Chebyshev Semi-iterative Method with Preconditioner**

For the system of linear equations **Ah** = **b** obtained by numerically discreting Richards' equation, the left preprocessing method is given by:

$$\mathbf{M}^{-1}\mathbf{A}\mathbf{h} = \mathbf{M}^{-1}\mathbf{b} \tag{4.1}$$

Consequently, the condition number of **M**−1**A** is greatly reduced.

Let **A** be a non-singular matrix. Matrix **A** can be split into **A** = **D** − **Q**, with **D** a non-singular matrix and is the diagonal matrix of matrix **A**. Then, if **H** = **D**−1**Q** has spectral radius less than 1.0, one obtains:

$$\mathbf{A}^{-1} = \left(\sum\_{k=0}^{\infty} \mathbf{H}^k\right) \mathbf{D}^{-1} \tag{4.2}$$

Matrix **M** is defined as:

$$\mathbf{M} = \mathbf{D} (\mathbf{I} + \mathbf{H} + \dots + \mathbf{H}^{m-1})^{-1} \tag{4.3}$$

where *m* is a natural number. Furthermore, matrix **M** is considered as a preconditioner to estimate the condition number of **M**−<sup>1</sup>**A**. From Eq. (4.3) and **A** = **D**−**Q**, one can obtain:

$$\mathbf{M}^{-1}\mathbf{A} = \left(\mathbf{I} + \mathbf{H} + \dots + \mathbf{H}^{m-1}\right)\mathbf{D}^{-1}(\mathbf{D} - \mathbf{Q}) = \mathbf{I} - \mathbf{H}^{m} \tag{4.4}$$

If λ1,...,λ*n* are the eigenvalue of matrix **H**, the eigenvalues of **M** are 1 <sup>−</sup> <sup>λ</sup>*<sup>m</sup> i* . If **D** and **A** are symmetric positive-definite matrices, then the eigenvalues of **H** are real. Assuming that λ1 ≤ ··· ≤ λ*n* < 1, then one can have:

$$\text{cond}(\mathbf{M}^{-1}\mathbf{A}) = \frac{\lambda\_{\text{max}}(\mathbf{M}^{-1}\mathbf{A})}{\lambda\_{\text{min}}(\mathbf{M}^{-1}\mathbf{A})} = \begin{cases} \frac{1-\lambda\_{1}^{m}}{1-\lambda\_{n}^{m}}, \text{ if } \lambda\_{1} \ge 0 \text{ or } \lambda\_{1} < 0, \ m \text{ is odd},\\ \frac{1-\sigma^{m}}{1-\lambda\_{n}^{m}}, \text{ if } \lambda\_{1} < 0, \ |\lambda\_{n}| \ge |\lambda\_{1}|, \ m \text{ is even},\\ \frac{1-\sigma^{m}}{1-\lambda\_{1}^{m}}, \text{ if } \lambda\_{1} < 0, \ |\lambda\_{1}| > |\lambda\_{n}|, \ m \text{ is even}, \end{cases} \tag{4.5}$$

where σ = min*i*|λ*i*|.

From Eq. (4.5), the condition number of matrix **M**−1**A** decreases with increasing *m*. Hence, the preprocessing can effectively decrease the condition number of this matrix and enhance the convergence rate.

Based on the Chebyshev polynomial acceleration process and the Gauss– Seidel iterative method of preconditioning, this chapter develops an improved preconditioning Chebyshev semi-iterative method (P-CSIM). First, the improved Gauss–Seidel iterative format using the preconditioning (Eq. 4.3) is written as:

$$\mathbf{h}^{k+1} = \mathbf{G}\_{\mathbf{M}} \mathbf{h}^{k} + \mathbf{b}\_{\mathbf{M}} \tag{4.6}$$

where **G**M denotes an iterative matrix and **s**M denotes a column matrix of preprocessing. Furthermore, the error vector for the general polynomial acceleration process of Eq. (4.6) is expressed as:

$$
\varepsilon^k = \mathcal{Q}\_k(\mathbf{G\_M}) \varepsilon^0 \tag{4.7}
$$

where *Qk* (**G**M) <sup>≡</sup> *ak*,0**I**+*ak*,1**G**<sup>M</sup> +···+*ak*,*k***G***<sup>k</sup>* M is a matrix polynomial, and the only additional condition being Σ*<sup>k</sup> <sup>i</sup>*=<sup>0</sup>*ak*,*i* = 1. The virtual spectral radius of the matrix *Qk* (**G**M) is then defined as:

$$\overline{\mathbf{S}}(\mathcal{Q}\_k(\mathbf{G}\_\mathcal{M})) = \max\_{m(\mathbf{G}\_\mathcal{M}) \le x \le M(\mathbf{G}\_\mathcal{M})} |\mathcal{Q}\_k(x)| \tag{4.8}$$

where *M*(**G**M) and *m*(**G**M) denote the algebraic maximum and minimum eigenvalues of **G**M, respectively.

It can be proved that the matrix polynomial *Qk* (**G**M) that ensures **S**(*Qk* (**G**M)) reaches a minimum is unique and can be defined by a Chebyshev polynomial. Furthermore, it is needed to find the polynomial *Pk* (*x*), for which *Pk* (1) = 1, and satisfies the following conditions (Theorem 4.2.1 of Hagemam and Young 1981):

$$\max\_{m(\mathbf{G}\_{\mathbf{M}}) \le x \le M(\mathbf{G}\_{\mathbf{M}})} |P\_k(\mathbf{x})| \le \max\_{m(\mathbf{G}\_{\mathbf{M}}) \le x \le M(\mathbf{G}\_{\mathbf{M}})} |Q\_k(\mathbf{x})| \tag{4.9}$$

where *Qk* (*x*) satisfies *Qk* (1) = 1. Assume that the eigenvalues of **G**M are real numbers, *m*(**G**M) = α, and *M*(**G**M) = β. According to Theorem 2.2.1 of Hagemam and Young (1981, p. 21), let:

$$
\alpha = \lambda\_1 \le \lambda\_2 \le \cdots \le \lambda\_n = \beta < 1, \quad \alpha \ne \beta \tag{4.10}
$$

### 4.2 A Chebyshev Semi-iterative Method with Preconditioner 75

To obtain *Pk* (*x*), a linear transformation is required as follows:

$$r = \frac{2x - \alpha - \beta}{\beta - \alpha} \tag{4.11}$$

Moreover, one can obtain:

$$\begin{cases} \left. r \right|\_{x=\alpha} = -1, \left. r \right|\_{x=\beta} = 1, \\ \left. x = \frac{(\beta - \alpha)r + (\beta + \alpha)}{2} \right. \\ \left. \mathcal{Q}\_k(\mathbf{x}) = \mathcal{Q}\_k \left( \frac{(\beta - \alpha)r + (\beta + \alpha)}{2} \right) =: P\_k(r). \end{cases} \tag{4.12}$$

A new parameter *z*<sup>∗</sup> is introduced as follows:

$$z^\* = r|\_{x=1} = \frac{2 - \alpha - \beta}{\beta - \alpha} \tag{4.13}$$

Note that *z*<sup>∗</sup> > 1 and *Pk* (*z*∗) = *Qk* (1) = 1. Here, *Pk* (*x*) can be defined as:

$$P\_k(\mathbf{x}) \equiv T\_k(r) / T\_k(z^\*) \tag{4.14}$$

Additionally, the *k*th Chebyshev polynomials on *Tk* (*r* ) for any non-negative integer *k* and *r* can be defined by the following 3-term recurrence equation (Arioli and Scott 2014):

$$\begin{cases} T\_0(r) = 1, & T\_1(r) = r, \\ T\_{k+1}(r) = 2rT\_k(r) - T\_{k-1}(r), & k \ge 1. \end{cases} \tag{4.15}$$

Using Theorem 4.2.1 (Hagemam and Young 1981), it is easy to prove that polynomials *Pk* (*x*) satisfy the recurrence relation. Therefore, according to Eqs. (4.6), (4.14), and (4.15) and Theorem 3.2.1 (Hagemam and Young 1981, p. 41), the 3-term recurrence equations for the P-CSIM are written as:

$$\begin{cases} \mathbf{h}^{1} = \overline{\boldsymbol{\gamma}} (\mathbf{G}\_{\mathbf{M}} \mathbf{h}^{0} + \mathbf{b}\_{\mathbf{M}}) + (1 - \overline{\boldsymbol{\gamma}}) \mathbf{h}^{0}, \\\mathbf{h}^{k+1} = \overline{\boldsymbol{\rho}}\_{k+1} [\overline{\boldsymbol{\gamma}} (\mathbf{G}\_{\mathbf{M}} \mathbf{h}^{k} + \mathbf{b}\_{\mathbf{M}}) + (1 - \overline{\boldsymbol{\gamma}}) \mathbf{h}^{k}] + (1 - \overline{\boldsymbol{\rho}}\_{k+1}) \mathbf{h}^{k-1} k \ge 1, \end{cases} \tag{4.16}$$

in which

$$\overline{\gamma} = \frac{2}{2 - \beta - \alpha} \tag{4.17}$$

$$\overline{\rho}\_{k+1} = \frac{2r(1)T\_k(r(1))}{T\_{k+1}(r(1))}\tag{4.18}$$

Equation (4.16) can be further expressed as:

76 4 Improved Linear and Nonlinear Iterative Methods for Rainfall …

$$\begin{cases} \mathbf{h}^{1} = \frac{2}{2 - \beta - \alpha} \left( \mathbf{G}\_{\mathbf{M}} \mathbf{h}^{0} + \mathbf{b}\_{\mathbf{M}} \right) - \frac{a + \beta}{2 - \beta - \alpha} \mathbf{h}^{0}, \\\ \mathbf{h}^{k+1} = 2 \frac{2 \mathbf{G}\_{\mathbf{M}} - (a + \beta) \mathbf{I}}{\beta - \alpha} \frac{T\_{k}(z^{\*})}{T\_{k+1}(z^{\*})} \mathbf{h}^{k} - \frac{T\_{k-1}(z^{\*})}{T\_{k+1}(z^{\*})} \mathbf{h}^{k-1} + \frac{4 \mathbf{b}\_{\mathbf{M}}}{\beta - \alpha} \frac{T\_{k}(z^{\*})}{T\_{k+1}(z^{\*})} \quad k \ge 1, \end{cases} \tag{4.19}$$

where *Tk* (*z*∗) is recursively calculated using Eq. (4.15); here, **I** is an identity matrix.

From Eq. (4.19), the convergence rate of P-CSIM mainly depends on the accuracy in estimating the maximum and minimum eigenvalues β and α of the iterative matrix **G**M. This chapter assumes 0 < −α = β < 1 based on the basic assumptions of Hagemam and Young (1981). By the above processing, the proposed P-CSIM achieves a better convergence and computational efficiency in solving the linear algebraic equations.

A brief implementation flowchart of P-CSIM is described in Fig. 4.1 (Zhu et al. 2020). To compare the speed-up of the improved method, the speed-up ratio is defined as:

$$\mathcal{S}\_{\text{GS}/\text{P}} = \frac{T\_{\text{GS}}}{T\_{\text{P-CSIM}}} \tag{4.20}$$

where *T*GSIM and *T*P-CSIM represent the computing times for GSIM and P-CSIM.

**Fig. 4.1** Brief flowchart of P-CSIM

### *4.2.1 Numerical Tests*

The mathematical model, boundary conditions, and unsaturated soil parameters of the numerical tests are consistent with Chap. 3, and the analytical solutions are expressed as Eqs. (2.3)–(2.5). The total duration of the simulation was 5 h, and the convergence criterion was set to 10−8 for the different iteration schemes. In preprocessing, the expansion order *m* is set to 20. To verify the accuracy, efficiency, and robustness of the proposed method, the time steps were set to 0.01 h, 0.005 h, 0.002 h, and 0.001 h, respectively, and the grid sizes were set to 0.4 m, 0.2 m, and 0.1 m.

Figure 4.2 demonstrates the relationship between expansion order *m* and the condition number for different time steps in this test. The condition number of coefficient matrix *A* increases with the time step. However, the condition number can be decreased to 1.0 through the preprocessing method. Additionally, Fig. 4.2a– d indicates that the condition number is close to 1.0 when *m* = 6, indicating that preprocessing effectively improves the ill-conditioning of the matrix.

Table 4.1 sums up the iterations, the computing time obtained from the different iterative schemes with different grid sizes, and the time steps in this test. Compared with other methods, the improved P-CSIM has the lowest iterations and the least computational cost at different time steps and grid sizes. When Δ*z* = 0.1, the

**Fig. 4.2** Relationship between parameter *m* and condition number at different time steps

iterations of the P-CSIM are much smaller than that of other methods (Fig. 4.3a). Moreover, the speed-up ratio of P-CSIM relative to the other three methods at Δ*t* = 0.005 (Fig. 4.3b) shows that it is the largest with JIM, CSIM, and GSIM following in sequence. This result also demonstrates that the proposed P-CSIM represents a vast improvement over CSIM.

RMSE reaches an order of 10−5 at *t* = 5 h (Fig. 4.4) and decreases with decreasing grid size for different iterative methods. Moreover, the accuracy of the proposed P-CSIM shows a slight improvement compared with that of GSIM and CSIM.

The results illustrate that the improved P-CSIM has better computational efficiency, stability, and accuracy than the other three iterative methods.


**Table 4.1** Numerical results of different iterative methods (Zhu et al. 2020)

**Fig. 4.3** Numerical results of different iterative methods at *t* = 0.005 h: **a** iterations and **b** speed-up ratio

**Fig. 4.4** Root mean squared error (RMSE) associated with several grid sizes at *t* = 5 h for the different iterative methods

### **4.3 Improved Gauss–Seidel Method**

When the first-order linear stationary iterative methods (Eq. 4.6) or other iterative schemes to solve linear equations are used, the iterative process may or may not converge. Even if it does, it may converge very slowly. In either condition, this chapter hopes to find an improved method to make non-convergent formats converge, and slow-convergent formats converge faster. The commonly used acceleration methods include extrapolation method, integral correction method, Chebyshev polynomial acceleration method, and preconditioning method.

### *4.3.1 Integral Correction Method*

The integral correction method is similar to Anderson acceleration method (Walker and Ni 2011; Both et al. 2018). Firstly, the mutually different approximate solutions **h**1, **h**2,… **h***c* (*c* > 1) can be obtained according to the basic iterative method, and **Ah***<sup>i</sup>* /= **b** (*i* = 1, 2, …, *c*). Then the information provided by **h***i* can be used to construct the vector **h** so that it is closer to the exact solution of the linear equations than **h***i* , that is, there is:

$$\|\|\mathbf{b} - \mathbf{A}\mathbf{h}\|\|\_{2} < \min\_{1 \le i \le c} \|\|\mathbf{b} - \mathbf{A}\mathbf{h}\_{i}\|\|\_{2} \tag{4.21}$$

The vector **h** that satisfies Eq. (4.21) is the correction solution of the linear equations with respect to **h**1, **h**2,… **h***c*. The process of determining **h** is the correction process, then the integral correction model can be expressed as:

$$\begin{cases} \mathbf{h} = \sum\_{i=1}^{c} \beta\_i \mathbf{h}\_i, & \sum\_{i=1}^{c} \beta\_i = 1, \\\ \|\mathbf{b} - \mathbf{A}\mathbf{h}\|\_2 = \min \end{cases} \tag{4.22}$$

Select arbitrarily *s* ∈ {1, 2,..., *c*}, then

$$\begin{aligned} \mathbf{h} &= \mathbf{h}\_s + \sum\_{i \neq s} \beta\_i (\mathbf{h}\_i - \mathbf{h}\_s), \\\\ \mathbf{r}(\mathbf{h}) &= \mathbf{b} - \mathbf{A} \mathbf{h} = (\mathbf{b} - \mathbf{A} \mathbf{h}\_s) - \sum\_{i \neq s} \beta\_i (\mathbf{A} \mathbf{h}\_i - \mathbf{A} \mathbf{h}\_s) = \mathbf{r}\_s + \sum\_{i \neq s} \beta\_i (\mathbf{r}\_i - \mathbf{r}\_s) \end{aligned} \tag{4.23}$$

Denote

$$\begin{aligned} \mathfrak{G}\_i &= \mathbf{r}\_i - \mathbf{r}\_s, \quad \mathbf{Q}\_s = \begin{bmatrix} \mathfrak{G}\_1, \dots, \mathfrak{G}\_{s-1}, \mathfrak{G}\_{s+1}, \dots, \mathfrak{G}\_c \end{bmatrix}, \\ \mathbf{y}\_s &= (\beta\_1, \dots, \beta\_{s-1}, \beta\_{s+1}, \dots, \beta\_c)^\top \end{aligned}$$

Therefore, Eq. (4.23) can be rewritten as **r**(**h**) = **r***s* + **Q***s***y***s*, so that,

$$\mathbf{h} = \sum\_{i=1}^{c} \beta\_i \mathbf{h}\_i = \beta\_1 \mathbf{h}\_1 + \beta\_2 \mathbf{h}\_2 + \dots + \beta\_c \mathbf{h}\_c \tag{4.24}$$

makes <sup>|</sup>**r**(**h**)|2 <sup>=</sup> min, which is equivalent to **y***s* <sup>∈</sup> <sup>C</sup>*c*−1 makes <sup>|</sup>**r***s* <sup>+</sup> **<sup>Q</sup>***s***y***s*|2 <sup>=</sup> min, and its minimal norm solution can be expressed as:

$$\mathbf{y}\_s = -\mathbf{Q}\_s^\dagger \mathbf{r}\_s \tag{4.25}$$

where **Q**† *<sup>s</sup>*is the generalized inverse matrix of **Q***s*. The parameter β*s* can be obtained by β*s* = 1 − Σ *<sup>i</sup>*/=*<sup>s</sup>*β*i* , and then the correction solution **h** can be determined using Eq. (4.24). In particular, for *c* = 2, one can gain:

$$\begin{cases} \beta\_1 = \mathbf{y}\_2 = -\mathbf{Q}\_2^\dagger \mathbf{r}\_2 = -\boldsymbol{\delta}\_1^\dagger \mathbf{r}\_2 = -\frac{(\mathbf{r}\_1 - \mathbf{r}\_2)^\prime \mathbf{r}\_2}{\|\mathbf{r}\_1 - \mathbf{r}\_2\|\_2^2} \\ \beta\_2 = 1 - \beta\_1 = \frac{(\mathbf{r}\_1 - \mathbf{r}\_2)^\prime \mathbf{r}\_1}{\|\mathbf{r}\_1 - \mathbf{r}\_2\|\_2^2} \end{cases} (\mathbf{s} = 2) \tag{4.26}$$

$$\begin{cases} \beta\_2 = \mathbf{y}\_1 = -\mathbf{Q}\_1^\dagger \mathbf{r}\_1 = -\boldsymbol{\delta}\_2^\dagger \mathbf{r}\_1 = -\frac{(\mathbf{r}\_2 - \mathbf{r}\_1)^H \mathbf{r}\_1}{\|\mathbf{r}\_2 - \mathbf{r}\_1\|\_2^2} \\ \beta\_1 = 1 - \beta\_2 = \frac{(\mathbf{r}\_2 - \mathbf{r}\_1)^H \mathbf{r}\_2}{\|\mathbf{r}\_2 - \mathbf{r}\_1\|\_2^2} \end{cases} (\mathbf{s} = 1) \tag{4.27}$$

It can be seen that the results of Eqs. (4.26) and (4.27) are consistent. Therefore, the improved Gauss–Seidel iterative method (IC(*c*)-GS) based on the integral correction method can be summarized as follows:


### *4.3.2 Multistep Preconditioner Method*

For a given system of linear equations, **Ah**<sup>∗</sup> = **b**, the condition number may be used to evaluate whether a given non-singular matrix **A** is ill-conditioned. The condition number for a non-singular matrix **A** can be defined as:

$$\text{Cond}(\mathbf{A}) = \|\mathbf{A}\| \cdot \|\mathbf{A}\|^{-1} \tag{4.28}$$

where the matrix norm is the Frobenius norm.

The preconditioning method can be an acceleration technique for the iterative solution of the system of linear equations. It can greatly improve the ill-condition of the original system of linear equations, thereby accelerating the convergence rate of the iterative method (Benzi 2002). The preconditioning process is usually to find the matrix **M** (preconditioner) for the linear equations **Ah**<sup>∗</sup> = **b**. There are many choices and forms of **M**. The preconditioning on the left for **Ah**<sup>∗</sup> = **b** can be written as:

$$\mathbf{M}^{-1}\mathbf{A}\mathbf{h}^\* = \mathbf{M}^{-1}\mathbf{b} \tag{4.29}$$

**Ah**<sup>∗</sup> = **b** is preconditioned from the right:

$$\mathbf{A}\mathbf{M}^{-1}\mathbf{y} = \mathbf{b}, \quad \mathbf{h}^\* = \mathbf{M}^{-1}\mathbf{y} \tag{4.30}$$

The preconditioning on both sides of **Ah**<sup>∗</sup> = **b** is:

$$\mathbf{M}\_1^{-1}\mathbf{A}\mathbf{M}\_2^{-1}\mathbf{z} = \mathbf{M}\_1^{-1}\mathbf{b}, \quad \mathbf{h}^\* = \mathbf{M}\_2^{-1}\mathbf{z}, \quad \mathbf{M} = \mathbf{M}\_1\mathbf{M}\_2\tag{4.31}$$

Based on the preconditioner increasing the complexity of the iterative method, here the preconditioner should be simple to construct and greatly improve the convergence rate of the iterative method. Therefore, the left preconditioning (Eq. 4.29) is adopted in this chapter. A symmetric Gauss–Seidel (SGS) preconditioner is written as:

### 82 4 Improved Linear and Nonlinear Iterative Methods for Rainfall …

$$\mathbf{M}\_{\rm SGS} = (\mathbf{D} - \mathbf{L})\mathbf{D}^{-1}(\mathbf{D} - \mathbf{U})\tag{4.32}$$

It is easy to know that combining the split matrix of the GS can quickly construct the preconditioner (Eq. 4.32). This preconditioning technique has almost no construction cost and is easy to implement. However, the result of one-step preconditioning is sometimes not good enough. Furthermore, a multistep SGS preconditioner is proposed to further obtain a sufficiently good and easier-to-solve system of the linear equations:

$$\prod\_{i=1}^{m} \mathbf{M}\_{\text{SGS}}^{-1}{}^{i} \mathbf{A} \mathbf{h}^{\*} = \prod\_{i=1}^{m} \mathbf{M}\_{\text{SGS}}^{-1}{}^{i} \mathbf{b} \tag{4.33}$$

where *m* represents the number of steps of the multistep preconditioner. Generally, a multistep preconditioner with a small number of steps that can lead to good convergence results in the iterative method, and a multistep preconditioner with a large number of steps may be more robust, but it may also cause higher computational cost of the iterative method. The empirical number of steps should be 2, 3, or 4, which can provide the preconditioning iterative method with good convergence rate. The Gauss–Seidel iterative method with multistep preconditioner (MP(*m*)-GS) can be summarized as follows:


For (*i* = 1; *i* ≤ *m*; *i*++). Compute an SGS preconditioner **M**−<sup>1</sup> SGS *<sup>i</sup>*using **A***i* . Drop the small entries in the matrix **M**−<sup>1</sup> SGS *<sup>i</sup>*relative to the parameter τ . Compute **A***i*+1 <sup>=</sup> **<sup>M</sup>**−<sup>1</sup> SGS *i* **A***i* . End the cycle.


The parameter τ can ensure the sparsity of the preconditioning matrix, thereby reducing the computational cost of the preconditioning process (Wang and Zhang 2003). In this chapter, the parameter τ is uniformly set to 10−4. It can be easily found that when *m* = 0, MP(*m*)-GS completely degenerates to GS.

### *4.3.3 Mixed Method*

The integral correction method utilizes the solution vector error correction to make the iterative method obtain a faster convergence rate, while the multistep preconditioning method preprocesses the original system of linear equations into a more easily solved system of linear equations, making the subsequent iterative method more efficient. The two methods are different. Therefore, the two methods are combined to propose an improved Gauss–Seidel iterative method (ICMP(*m*)-GS) with multistep preconditioner based on the integral correction method. The method is summarized as follows:


To verify the calculation efficiency, the speed-up ratio is defined as:

$$S\_{\rm GS/MP(m)\cdot\rm GS} = \frac{T\_{\rm GS}}{T\_{\rm MP(m)\cdot\rm GS}}\tag{4.34}$$

where *T*GS and *T*MP(*m*)-GS are the runtimes of GS and MP(*m*)-GS, respectively.

### *4.3.4 Numerical Tests*

### **4.3.4.1 1D Steady-State Unsaturated Infiltration**

The mathematical model is shown in Fig. 1.1. The model parameters are: α = 8 × 10−3, θs = 0.35, θr = 0.14, and *h*d = − 103 m. For the integral correction method, to test the influence of the parameter *c* on the improved iterative method, *c* is set to 2, 5, 10, and 20. For the multistep preconditioning method, *m* is 1, 2, 3, and 4 to verify the effectiveness of the algorithm. In Fig. 4.5, the change in the space step Δ*z* has a small impact on the acceleration effect of the IC(*c*)-GS iterative method, while the change of the parameter *c* has a greater impact on it. When the parameter *c* is less than 10, the speed-up ratio of IC(*c*)-GS relative to the GS is relatively small. And when *c* ≥ 10, the speed-up ratio increases greatly. It can be found that the acceleration effects of *c*  = 10 and *c* = 20 are roughly similar. Therefore, the parameter *c* of the IC(*c*)-GS is set to 10 in this chapter. In Fig. 4.6, the convergence rate of the Jacobi is the slowest,

**Fig. 4.5** Comparison of the speed-up ratio of IC(*c*)-GS relative to GS under different parameter *c*

followed by the GS, and the SOR is faster than the Jacobi and GS. It is easy to obtain that when *c* ≥ 10, the convergence rate of IC(*c*)-GS is faster than SOR. When *m* = 3 and 4, the convergence rate of MP(*m*)-GS is also faster than SOR. For the proposed mixed method ICMP(*m*)-GS, it can be found from the enlarged convergence graph that its convergence rate is faster than SOR, IC(*c*)-GS, and MP(*m*)-GS, regardless of the value of *m*, and increases with the increase of *m*.

In Fig. 4.7, the speed-up ratio of MP(*m*)-GS and ICMP(*m*)-GS relative to GS is very large when Δ*z* = 0.025 m, and the acceleration effect is much better than IC(*c*)- GS. The acceleration effect of ICMP(*m*)-GS is better than MP(*m*)-GS. Additionally, the speed-up ratio of ICMP(*m*)-GS relative to GS increases slowly when the parameter *m* is 3 and 4, which is related to the increase of the convergence rate. In other words, the improved method ICMP(*m*)-GS increases the computational complexity and cost as increasing *m*, but the number of iterations does not greatly decrease. The numerical results are listed in Table 4.2. The number of iterations and running time of GS and ICMP(3)-GS increase as Δ*z* decreases, but the number of iterations of ICMP(3)-GS is much smaller than that of GS, and the difference is about 10,000 times. The running time of ICMP(3)-GS is much shorter than that of GS, and the difference is about 100 times.

### **4.3.4.2 1D Transient Unsaturated Infiltration**

The boundary conditions can be written as follows:

$$h(z=0) = h\_{\mathbb{d}}\tag{4.35}$$

**Fig. 4.6** Comparison of convergence rates for the different iterative methods

$$h(z=10) = 0\tag{4.36}$$

The parameters are set to: α = 1 × 10−4, θs = 0.50, θr = 0.11, and the saturated permeability coefficient *k*s = 2.5 × 10−8 m/s. The total simulation time is 10 h. To verify the calculation efficiency of the proposed method, the space step is taken as 0.05 m, 0.025 m, and 0.0125 m, and the time step is 0.2 h, 0.1 h, and 0.05 h, respectively.

Figure 4.8 shows the changes in the condition number of the coefficient matrix **A**  with and without the algorithm MP(*m*) under different numerical discrete conditions. The condition number of the coefficient matrix increases with the decrease of the

**Fig. 4.7** Effect of parameter *m* on the improved iterative methods MP(*m*)-GS and ICMP(*m*)-GS for Δ*z* = 0.025 m


**Table 4.2** Numerical solutions of steady-state infiltration

space step and decreases with the decrease of the time step. When the preconditioning method is not used, the condition number reaches in the order of 104. However, the condition number can be greatly reduced using the algorithm MP(*m*) and decrease with the increase of the parameter *m*. Particularly after using algorithm MP(4), the condition number even approaches 1.0, which indicates that the multistep preconditioning process can effectively improve the ill-condition of the system of the linear equations.

Tables 4.3 and 4.4 list the average number of iterations and running time of the conventional iterative methods GS and SOR and the improved method ICMP(*m*)- GS. First of all, it can be seen that the performance of SOR is better than that of GS. However, the number of iterations of SOR is not stable under a smaller space step. Compared with GS and SOR, the number of iterations and running time of the improved ICMP(*m*)-GS is less than that of the conventional method. The number of iterations increases with the decrease of Δ*z* and decreases with the decrease of Δ*t*. The running time increases as Δ*z* and Δ*t* decrease. At the same time, the

**Fig. 4.8** Condition number with and without the MP algorithm for different conditions **a** Δ*z* = 0.05 m and **b** Δ*z* = 0.025 m


**Table 4.3** Number of iterations for solving steady-state infiltration

performance of ICMP(*m*)-GS with *m* = 3 and 4 is much better than that of GS and SOR, and it has a very significant acceleration effect. Besides, the acceleration effect of the improved method ICMP(*m*)-GS also has certain rules under different grid sizes. In Fig. 4.9, with the decrease of the space step Δ*z* and the increase of the time step Δ*t*, the speed-up ratio of ICMP(3)-GS relative to GS has a tendency to increase. In Fig. 4.10, the transient numerical solutions obtained by the ICMP(*m*)-GS method are also very consistent with the analytical solutions. This result further verifies that the proposed method has a faster convergence rate and a higher acceleration effect than the conventional methods.

### **4.3.4.3 2D Transient Unsaturated Infiltration**

The geometry and boundary conditions for this example are illustrated in Fig. 2.2, where *L* = 1 m and *W* = 1 m. The soil was assumed to be silt. The parameters of soil were described by Liu et al. (2015) and included θs = 0.35, θr = 0.14, α = 8 × 10−3, and *k*s = 9 × 10−4 m/h. The boundary conditions are as follows:


**Table 4.4** Running time for solving steady-state infiltration

**Fig. 4.9** Comparison of the speed-up ratios of ICMP(3)-GS relative to GS under different grid sizes

$$h(0, z, t) = h\_{\mathbb{d}} \tag{4.37}$$

$$h(W, z, t) = h\_{\mathbb{d}} \tag{4.38}$$

$$h(\mathbf{x}, \mathbf{0}, t) = h\_{\mathbf{d}} \tag{4.39}$$

$$h(\mathbf{x}, L, t) = \ln\left( \left( 1 - \mathbf{e}^{ah\_d} \right) \sin\left(\frac{\pi \chi}{W}\right) + \mathbf{e}^{ah\_d} \right) / \alpha \tag{4.40}$$

**Fig. 4.10** Comparison of the numerical solutions obtained using ICMP(*m*)-GS with analytical solutions

The total simulation time and time step are set to 5 h and 0.1 h, respectively. The pressure profile computed using ICMP(*m*)-GS at *t* = 5 h is presented in Fig. 4.11, which is consistent with the analytical solution of this issue.

In Fig. 4.12, the condition number of the coefficient matrix can be greatly reduced using the proposed algorithm MP(*m*) and decreases with the increase of *m*. Figure 4.13a demonstrates the number of iterations of ICMP(*m*)-GS and GS. The number of iterations of ICMP(*m*)-GS is much smaller than that of the GS and decreases as *m* decreases. For *m* > 2, the number of iterations decreases slowly. Figure 4.13b indicates that the running time of ICMP(*m*)-GS is less than that of the GS, and the difference between the running time of GS and ICMP(3)-GS is about 10 times. Additionally, it can be found that when Δ*x*, Δ*z* = 0.025 m, the running time of ICMP(3)-GS is slightly shorter than that of ICMP(4)-GS. This is because when *m* = 3, the condition number of the coefficient matrix has been greatly reduced (Fig. 4.12). Thus increasing *m* cannot improve the computation efficiency. This result further demonstrates that the proposed method ICMP(*m*)-GS can achieve a significant acceleration compared to the GS, particularly at a smaller grid size.

**Fig. 4.11** Numerical solutions obtained using ICMP(*m*)-GS for 2D transient infiltration problem at 5 h

**Fig. 4.12** Relationship between the condition number and the parameter *m* for different grid sizes

**Fig. 4.13** Numerical results for this test: **a** number of iterations and **b** running time

### **4.4 Nonlinear Iterative Methods and Improvements**

### *4.4.1 Newton and Picard Methods*

In Fig. 4.14, the Richards' equation needs to be modified to satisfy rainfall infiltration into unsaturated soil slopes (Iverson 2000):

$$\frac{\partial}{\partial z}\left[K(h)\left(\frac{\partial h}{\partial z} + \cos \beta\right)\right] = \frac{\partial \theta}{\partial t} \tag{4.41}$$

where β is the slope angle.

**Fig. 4.14** Schematic map of control volume method based on non-uniform grid nodes and 1D infiltration model for homogeneous soil slope

Finite volume method is used to discretize Eq. (4.41). First, the interval on the *z*-axis is divided into *N* equal parts, and the simulation time is divided into *M* equal parts. Furthermore, Eq. (4.41) is integrated as follows (Patankar 1980):

$$\int\_{i}^{i+\Delta z} \int\_{i}^{i+\Delta t} \frac{\partial \theta}{\partial t} \mathrm{d}t \,\mathrm{d}z = \int\_{i}^{i+\Delta t} \int\_{i}^{i} \frac{\partial}{\partial z} \left[ K(h) \left( \frac{\partial h}{\partial z} + \cos \beta \right) \right] \mathrm{d}z \,\mathrm{d}t \tag{4.42}$$

By simplification, one can gain:

$$\frac{\mathbf{K}\_{i+1/2}^{j}\left(h\_{i+1}^{j} - h\_{i}^{j}\right)}{\Delta z} - \frac{\mathbf{K}\_{i-1/2}^{j}\left(h\_{i}^{j} - h\_{i-1}^{j}\right)}{\Delta z} + \left(\mathbf{K}\_{i+1/2}^{j} - \mathbf{K}\_{i-1/2}^{j}\right)\cos\beta$$

$$= \Delta z \mathbf{C}\_{i}^{j-1/2} \frac{h\_{i}^{j} - h\_{i}^{j-1}}{\Delta t}, \quad 1 \le i \le N - 1, \ 1 \le j \le M + 1\tag{4.43}$$

where *i* represents discrete nodes along the *z*-axis (except boundary nodes) and *j*  represents time nodes. *C* is the specific moisture capacity, defined as *C*(*h*) = ∂θ/∂*h*; *Ki*+1/2 and *Ki*−1/2 represent the harmonic average of the hydraulic conductivity corresponding to adjacent nodes.

Newton and Picard iterative methods to solve Eq. (4.43) are a popular linearization technique. Newton method based on the one-dimensional Richards' equation in the form of pressure head (*h*) is expressed as:

$$\left(\mathbf{A}^{j,k} + \frac{\partial \mathbf{A}^{j,k}}{\partial \mathbf{h}^{j,k}} \mathbf{h}^{j,k} + \frac{\mathbf{C}^{j,k}}{\Delta t} + \frac{1}{\Delta t} \frac{\partial \mathbf{C}^{j,k}}{\partial \mathbf{h}^{j,k}} (\mathbf{h}^{j,k} - \mathbf{h}^{j-1}) - \frac{\partial \mathbf{F}^{j,k}}{\partial \mathbf{h}^{j,k}}\right)$$

$$\left(\mathbf{h}^{j,k+1} - \mathbf{h}^{j,k}\right) = \mathbf{F}^{j,k} - \mathbf{A}^{j,k} \mathbf{h}^{j,k} - \frac{\mathbf{C}^{j,k}}{\Delta t} \left(\mathbf{h}^{j,k} - \mathbf{h}^{j-1}\right) \tag{4.44}$$

where *k* is the number of iterations and **A** is a tridiagonal matrix, at node *i*, which is given by:

$$\mathbf{A}\_{i} = \left[ \frac{K\_{i-1/2}^{j,k}}{\Delta z}, -\frac{K\_{i-1/2}^{j,k}}{\Delta z} - \frac{K\_{i+1/2}^{j,k}}{\Delta z}, \frac{K\_{i+1/2}^{j,k}}{\Delta z} \right] \tag{4.45}$$

**C** and **F** can be written as follows at node *i*:

$$\mathbf{C}\_{i} = -\Delta z \mathbf{C}\_{i}^{j-1/2,k} \tag{4.46}$$

$$\mathbf{F}\_i = \left( K\_{i-1/2}^{j,k} - K\_{i+1/2}^{j,k} \right) \cos \beta \tag{4.47}$$

The modified Picard method (PI), also known as the fixed-point method (Schrefler and Zhan 1993), can be expressed as:

4.4 Nonlinear Iterative Methods and Improvements 93

$$\left(\mathbf{A}^{j,k} + \frac{\mathbf{C}^{j,k}}{\Delta t}\right) \left(\mathbf{h}^{j,k+1} - \mathbf{h}^{j,k}\right) = \mathbf{F}^{j,k} - \mathbf{A}^{j,k}\mathbf{h}^{j,k} - \frac{\mathbf{C}^{j,k}}{\Delta t} \left(\mathbf{h}^{j,k} - \mathbf{h}^{j-1}\right) \quad (4.48)$$

It can be seen from Eq. (4.48) that the modified Picard method is a simplification of the Newton method, ignoring the derivative term in the Jacobian matrix. In addition, Eq. (4.48) can be further simplified into the standard Picard format:

$$\left(\mathbf{A}^{j,k} + \frac{\mathbf{C}^{j,k}}{\Delta t}\right)\mathbf{h}^{j,k+1} = \mathbf{F}^{j,k} + \frac{\mathbf{C}^{j,k}}{\Delta t}\mathbf{h}^{j-1} \tag{4.49}$$

Newton's method is quadratic convergent (Li 1993), but needs to calculate the first derivatives of matrices **A** and **C** and vector **F**, and the modified Picard method only converges linearly (Li 1993). Due to Newton's method increases the algebraic complexity and computational cost of constructing derivative terms in Jacobian matrices, thus modified Picard method has been widely investigated for *h*-based and mixed form of Richards' equations. Furthermore, there is still a large computational cost due to the need to recalculate the iteration matrix for each iteration in the Picard method (Farthing and Ogden 2017). For some unfavorable numerical conditions, such as infiltration into dry soils and/or layered soils with very different permeability coefficients, this results in numerical difficulties or even non-convergence of conventional Picard methods (Zha et al. 2017).

Normally, the uniform discrete process does not have good numerical convergence under some unfavorable numerical conditions, such as rainfall infiltration into dry soils (Zha et al. 2017). Therefore, a non-uniform grid in the form of Chebyshev is adopted (Wu et al. 2020). The coordinates in the two-layer soils in Fig. 4.14 can be expressed as follows:

$$z\_{li} = \cos(i\pi/N\_1) \times \frac{L\_1}{2} + \frac{L\_1}{2}, \quad i = N\_1, N\_1 - 1, \dots, 0 \tag{4.50}$$

$$z\_{2i} = \cos(i\pi/N\_2) \times \frac{L\_2}{2} + \frac{L\_2}{2} + L\_1, \quad i = N\_2 - 1, N\_2 - 2, \dots, 0 \tag{4.51}$$

where *L*1 and *L*2 are the height of soil layers 1 and 2, respectively; *i* denotes nonuniform grid nodes; and *N*1 and *N*2 represent the number of nodes in layers 1 and 2, respectively. In addition, for the convenience of comparison and analysis, one can have *N* = *N*1+*N*2. The discretized equation is now derived by integrating Eq. (4.41) over the non-uniform finite volume in Fig. 4.14 and over the time interval from *t* to *t* + Δ*t*. Thus,

$$\int\_{i\_{\rm W}}^{i\_{\rm f}} \int\_{t}^{t+\Delta t} \frac{\partial \theta}{\partial t} \mathrm{d}t \,\mathrm{d}z = \int\_{t}^{t+\Delta t} \int\_{i\_{\rm W}}^{i\_{\rm F}} \frac{\partial}{\partial z} \left[ K(h) \left( \frac{\partial h}{\partial z} + \cos \beta \right) \right] \mathrm{d}z \,\mathrm{d}t \tag{4.52}$$

Furthermore, this chapter can obtain

$$\begin{split} \frac{K\_{i\mathbb{E}}^{j,k}\left(h\_{i+1}^{j,k+1} - h\_{i}^{j,k+1}\right)}{(\delta z)\_{\mathbb{E}}} - \frac{K\_{i\mathbb{w}}^{j,k}\left(h\_{i}^{j,k+1} - h\_{i-1}^{j,k+1}\right)}{(\delta z)\_{\mathbb{W}}} + \left(K\_{i\mathbb{E}}^{j,k} - K\_{i\mathbb{w}}^{j,k}\right)\cos\beta\\ = C\_{i}^{j-1/2,k}(\Delta z)\_{i} \frac{h\_{i}^{j,k+1} - h\_{i}^{j-1}}{\Delta t} \end{split}{\tag{4.53}}$$

According to Eqs. (4.43) and (4.53), the matrix form of the standard Picard iterative method can be abbreviated as follows:

$$\mathbf{A}(\mathbf{h}^{j,k})\mathbf{h}^{j,k+1} = \mathbf{f}(\mathbf{h}^{j,k}, \mathbf{h}^{j-1})\tag{4.54}$$

For the solution of Eq. (4.54), the system of linear equations is first derived and then solved by the basic linear solution methods. The Gaussian elimination method is usually effective. After solving Eq. (4.54) for the first time, the coefficient matrix in Eq. (4.54) is recalculated using this first solution, and the new linear equation is solved. The iteration process is terminated when the convergence measure is satisfied using *l*<sup>∞</sup> norm as:

$$\frac{\left\|\mathbf{h}^{j,k+1} - \mathbf{h}^{j,k}\right\|\_{\infty}}{\left\|\mathbf{h}^{j,k}\right\|\_{\infty}} < \varepsilon \tag{4.55}$$

The standard Picard method according to Eq. (4.54) can be abbreviated as N-PI. Generally, the Picard method combined with non-uniform grids has a slower convergence rate and lower calculation efficiency. To improve the convergence rate of N-PI, an improved Picard method based on non-uniform two-grid correction scheme (NTG-PI) is proposed in this book. In addition, the improved Picard method by the adaptive relaxation method (NAR-PI) is adopted for the comparative study.

### *4.4.2 Improved Picard Method*

### **4.4.2.1 Adaptive Relaxed Picard Iterative Method (NAR-PI)**

The adaptive relaxation method can effectively improve the computational efficiency of the Picard method (Smedt et al. 2010). First, when *k* > 1, the adaptive relaxation method can be expressed in the following form:

$$\mathbf{h}^{k+1} \leftarrow \mathbf{h}^k + \lambda\_k (\mathbf{h}^{k+1} - \mathbf{h}^k); \quad \lambda\_k \in (0, 2) \tag{4.56}$$

where λ*k* is the relaxation coefficient of the *k*th iteration, and its size can be adjusted according to the generalized angle between the current iteration step increment Δ**h**  and the previous iteration step increment Δ**h**:

$$\delta = \arccos\left(\frac{\Delta \mathbf{h} \cdot \Delta \overline{\mathbf{h}}}{\|\Delta \mathbf{h}\| \|\Delta \overline{\mathbf{h}}\|}\right) \tag{4.57}$$

The angle δ in Eq. (4.57) represents the convergence trend of the numerical solution. An acute angle indicates better convergence, and an obtuse angle indicates oscillation. Therefore, when δ is an acute angle, λ*k* should increase; when δ is an obtuse angle, λ*k* should be appropriately reduced. After preliminary trial calculation, better calculation results can be obtained in the following forms (Smedt et al. 2010):

$$\lambda\_k = \begin{cases} \sqrt{2}\lambda\_{k-1}, & \delta < \pi/4 \\ \lambda\_{k-1}, & \pi/4 \le \delta < \pi/2 \\ \lambda\_{k-1}/\sqrt{2}, & \delta \ge \pi/2 \end{cases} \tag{4.58}$$

Therefore, the adaptive relaxed Picard iterative method based on the non-uniform grid (NAR-PI) can be summarized as follows:


### **4.4.2.2 Improved Picard Method Based on the Two-Grid Correction Scheme (NTG-PI)**

The algebraic multigrid method can improve the computational efficiency of nonlinear iterative methods (Wang and Schrefler 2003). In this chapter, the nonuniform two-grid correction scheme (NTG) is adopted to improve the computational efficiency of the classical Picard method. In Fig. 4.15, the non-uniform two-grid correction scheme is described.

First, the solution vectors of the fine and coarse grids are assumed to be **h**f and **h**c , respectively. In NTG, the Gauss–Seidel iterative method is usually used as the basic linear solution method instead of the Gaussian elimination method. The implementation steps of NTG-PI are expressed in detail as follows (Zhu et al. 2022a, b):


**Fig. 4.15** Schematic drawing of non-uniform two-grid correction scheme

(3) Solve the following residual equation for the coarse grid:

$$\mathbf{A}(\mathbf{h}^c)\mathbf{e}^c = \mathbf{r}^c \tag{4.59}$$

(4) Interpolate the coarse-grid error **e**c to the fine grid and correct the fine-grid approximation solution by:

$$\mathbf{h}^{\text{f}\mu} \leftarrow \mathbf{h}^{\text{f}\mu} + \mathbf{I} \cdot \mathbf{e}^{\text{c}} \tag{4.60}$$

### 4.4 Nonlinear Iterative Methods and Improvements 97

where **I** denotes the interpolation operator from coarse to fine grids.

(5) Relax μ times using Gauss–Seidel relaxation for the fine grid with corrected vector **h**fμ, where μ usually takes a small natural number. μ is set to 1 in this chapter.

If **h**f<sup>μ</sup> converges according to Eq. (4.55), then stop; otherwise, go to step 1. Steps 1 and 5 can eliminate the high-frequency components of the iteration error, while steps 2–4 can eliminate the remaining low-frequency smooth components. Through multiple cycles of the above steps, it is expected to eliminate the iteration error as soon as possible. And this makes the iterative method more efficient. In addition, it is needed to further define the interpolation (**I**) and restriction (**R**) operators. In Fig. 4.15, the restriction operator for nodes *i* can be expressed as:

$$u\_i = \frac{(\Delta z)\_1}{(\Delta z)\_i} u\_{2i-1} + \frac{(\Delta z)\_2}{(\Delta z)\_i} u\_{2i} + \frac{(\Delta z)\_3}{(\Delta z)\_i} u\_{2i+1} \tag{4.61}$$

The interpolation operator can be written as:

$$
\mu\_{2i} = \mu\_i \tag{4.62}
$$

$$u\_{2i+1} = \frac{(\Delta z)\_2}{(\Delta z)\_1 + (\Delta z)\_2} u\_i + \frac{(\Delta z)\_1}{(\Delta z)\_1 + (\Delta z)\_2} u\_{i+1} \tag{4.63}$$

$$
\mu\_{2i+2} = \mu\_{i+1} \tag{4.64}
$$

### *4.4.3 Numerical Tests—1D Transient Unsaturated Infiltration*

This test simulates the 1D transient infiltration in homogeneous soils to verify the effectiveness of the proposed schemes. The mathematical model is described in Fig. 4.14. The boundary conditions and the parameters are kept unchanged. The total simulation time is 10 h. The parameter μ is set to 1. Additionally, the number of nodes is taken as 100, 200, and 400, and the time steps are taken as 0.1 h, 0.05 h, and 0.01 h, respectively.

Figure 4.16a represents the maximum relative error (MRE) obtained by different methods under different time steps when the number of nodes *N* = 200. It can be seen that the range of MRE obtained by NTG-PI is between 0.35 and 4.6%. When *t*  is less than 6 h, the MRE decreases as the time step decreases. However, the range of MRE obtained by PI is between 2.7 and 72%, which increases with the increase of *t* and is much larger than that obtained by NTG-PI. Figure 4.16b shows the MRE obtained by different methods under different number of nodes when the time step Δ*t* = 0.01 h. It can be found that the MRE obtained by PI increases over time *t*. When

**Fig. 4.16** Comparison of the maximum relative error of different methods under different numerical conditions: **a** time step and **b** number of nodes *N* 

*t* is greater than 2 h, the error becomes larger and larger, while the MRE obtained by NTG-PI only ranges from 0.18 to 8%. In addition, the MRE of the two methods decreases as the number of nodes *N* increases.

In Fig. 4.17, the numerical solutions obtained by the two methods under Δ*t* = 0.01 h and *N* = 200 are compared with the exact solutions. It can be found that there is a great deviation between the numerical solution obtained by PI and the exact solution, particularly after the time is greater than 4 h (Fig. 4.17a). On the contrary, the numerical solution obtained by NTG-PI is very consistent with the exact solution, and there is no large relative error as the time increases (Fig. 4.17b). When *t* = 10 h, the RSE of NTG-PI decreases with the increase of *N* and the decrease of Δ*t* in Table 4.5, which seems to be better than that of PI. And the RE of NTG-PI is much smaller than that of PI, and the difference is about 100 times. In addition, Fig. 4.18 depicts the convergence rate of the proposed methods. It can be found that the convergence rate of the proposed method NTG-PI is faster than that of the conventional method N-PI. As the number of nodes *N* changes, the convergence rate of NTG-PI does not change much.

In Table 4.6, the speed-up ratio of the improved method NTG-PI relative to N-PI is much greater than that of NAR-PI relative to N-PI, with a difference of nearly 10 times. Since the parameter μ is only 1, the computational efficiency of the proposed NTG-PI is higher. This test indicates that the proposed NTG-PI can obtain higher numerical accuracy with a smaller number of nodes and has a faster convergence rate while reducing the computational cost.

**Fig. 4.17** Comparison of the numerical solutions obtained using different methods with analytical solutions: **a** PI and **b** NTG-PI


**Table 4.5** Numerical accuracy of transient infiltration in homogeneous soils at *t* = 10 h

### **4.5 Conclusions**

This chapter first discusses the classical linear stationary iterative methods, such as Jacobi iterative method, Gauss–Seidel iterative method, and SOR iterative method. Then, combined with the preprocessing technology and error correction method, a series of improved iterative methods is proposed. The proposed methods have better computational performance in the simulation of rainfall infiltration and can obtain good application results. The main conclusions are as follows:

(1) Chebyshev semi-iterative method with polynomial preconditioner (P-CSIM) is developed. Compared with conventional iterative methods such as Jacobi and Gauss–Seidel methods, P-CSIM has less iterations and calculation time in

**Fig. 4.18** Comparison of convergence rates for the proposed schemes



simulating rainfall. Numerical results indicate significant speed-up, at least 50% increase compared to GSIM, and higher computational efficiency compared to CSIM and JIM.


from the linearized Richards' equation, but also obtain higher numerical accuracy. Numerical results show that compared with the conventional methods GS and SOR and the improved methods IC(*c*)-GS and MP(*m*)-GS, ICMP(*m*)-GS can improve the convergence rate to a greater degree and has higher computational efficiency and calculation accuracy. This method is relatively simple, and it has a good application prospect for simulating unsaturated infiltration.

(4) In the conventional Picard method, the Gaussian elimination method is usually used to solve the system of linear equations. Although the solution obtained by this method is theoretically accurate, with the increase of the number of discrete nodes *N*, the calculation amount increases in the order of *N*3, and the computational efficiency significantly reduced. In the improved method NTG-PI, the Gauss–Seidel iterative method is used instead of the Gaussian elimination method as the basic method, and its solution speed is faster. Numerical results demonstrate that the numerical solution obtained by the NTG-PI is in good agreement with the analytical solution under extremely dry initial conditions. Compared with the conventional Picard method, NTG-PI can achieve higher numerical accuracy with fewer discrete nodes, while simulating rainfall infiltration with faster convergence.

### **References**


Hagemam LA, Young DM (1981) Applied iterative methods. Academic Press, New York


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Chapter 5 Slope Stability Analysis Based on Analytical and Numerical Solutions**

### **5.1 Introduction**

Infiltration into soil slopes is a fundamental concern in civil engineering. Rainfall infiltration leads to changes in pore-water pressure and reduces matric suction in soils, making it one of the main triggers of slope failure (Rahimi et al. 2010; Ali et al. 2014; Wu et al. 2020). Slope instabilities caused by water infiltration are called rainfall-induced landslides (Xu and Zhang 2010; Wu et al. 2020).

There are three methods that can be taken to investigate the effect of rainfall infiltration on pore-water pressure or pressure head profile and hence on unsaturated soil slope stabilities. The approaches include numerical simulation, field monitoring, and analytical analysis (Zhan et al. 2013). A number of numerical studies have been conducted to investigate the hydrodynamic behaviors of soil slopes due to rainfall infiltration and to estimate the influence on the slope stability (Ng and Shi 1998; Iverson 2000; Collins and Znidarcic, 2004; Kim et al. 2004; Zhang et al. 2005; Garcia et al. 2011; Ali et al. 2014). Factors affecting the soil slope stability due to rainwater infiltration comprise the rainfall characteristics, the saturated permeability coefficient, the geometry of the slope, and the boundary and initial soil moisture conditions (Ali et al. 2014). Laboratory and field experiments have been performed to investigate variations in matric suction due to rainfall infiltration to improve the understanding of the mechanism of rainfall-induced soil slope failures (Wu et al. 2015, 2017, 2020). Field monitoring is helpful in the study of the effects of rainfall infiltration on the slope stability, but it is a costly procedure (Zhan et al. 2013).

The analytical methods involve theoretical infiltration and approximate infiltration models (Wu et al. 2022). In a theoretical infiltration model, the partial differential equation describing water infiltration in soils (i.e., the Richards' equation) is proposed based on continuum mechanics, and the solution can be obtained through integral transformation or numerical methods (Srivastava and Yeh 1991; Zhu et al. 2020, 2022).

The analytical method requires that certain assumptions be made regarding the closed-form equations that are derived. If the assumptions can be shown to be reasonable, then the analytical methodology becomes a simple and practical tool for studying possible pressure head distribution and the factor of safety of the slope.

Coupled poromechanical approaches incorporating the unsaturated soils have been employed to analytically and numerically examine the hydrodynamic behaviors of partially or completely saturated slopes due to rainfall (Chen et al. 2005). These coupled poromechanical studies are in view of various constitutive models, which represent the mechanical response of the unsaturated soils using nonlinear elastic stress–strain relationships (Cho and Lee 2001; Zhang et al. 2005) or elastoviscoplastic and elastoplastic models. The coupled hydromechanical approach has been combined with a slope stability analysis according to finite element methods to examine the response of an unsaturated soil slope subjected to a rainstorm and to consider the variability in the soil properties (Zhang et al. 2005). Consequently, using the improved method proposed in the previous chapters, the infiltration equation of the unsaturated soil slope is solved, and the numerical and analytical solutions are applied to the slope model to study the effect of rainfall infiltration on the slope stability.

The practical infiltration process is very complicated and affected by many factors such as soil profile and rainfall conditions and becomes difficult to be described accurately with theoretical formulations (D'Aniello et al. 2019; Srivastava et al. 2020).

Due to the nonlinearity of differential equations, numerical software packages are frequently needed to solve the complex problem of infiltration for predicting moisture movement and pore-water pressure change in unsaturated soils (Masoudian et al. 2019; Yang et al. 2018; Zhang et al. 2016).

Numerical approaches often suffer from convergence and mass balance problems and also are inefficient and expensive in many conditions. To provide simplified solutions to the complex infiltration issue, a number of theoretical models based on the wetting front have been proposed (Iverson 2000; Conte and Troncone, 2008). A series of analytical approaches (Srivastava and Yeh 1991; Iverson 2000; Wu et al. 2022) was developed to calculate the pressure head change during rainfall. Unfortunately, these models can only be applied under the assumption that the rainfall intensity is a constant. However, in practice, there are often cases where the rainfall intensity is a variable function depending on duration.

The objective of this chapter is to study the application of improved numerical methods to the slope stability assessment and to discuss the influence of various factors on the slope stability. Combined with the monitoring data of slope pore-water pressure in the Tung Chung area of Hong Kong, the book carried out a comparative study to verify the accuracy and practicability of the proposed improved method. The Xiaoba landslide in Guizhou, China, a typical rainfall-induced landslide, was investigated using numerical methods.

### **5.2 Application of Unsaturated Soil Slope Stability Under Rainfall**

### *5.2.1 Application of Chebyshev Spectral Method in Slope Stability Analysis*

For unsaturated soil slopes (Fig. 5.1), the governing equations are modified as follows:

$$\frac{\partial}{\partial z} \left[ K\_z(h) \left( \frac{\partial h}{\partial z} + \cos \beta \right) \right] = \frac{\partial \theta}{\partial t} \tag{5.1}$$

where β is the slope angle.

According to modified Eq. (2.34), the governing equation considering the hydromechanical coupling in soil slopes can be obtained:

$$\frac{\partial}{\partial z^\*} \left[ K(h) \frac{\partial h}{\partial z^\*} \right] + \frac{\partial K(h)}{\partial z^\*} \cos \beta = C(h) \frac{\partial h}{\partial t} + \frac{\alpha\_\text{c} \wp\_\text{w} (\theta - \theta\_\text{r})}{(\theta\_\text{s} - \theta\_\text{t}) F} \frac{\partial h}{\partial t} \tag{5.2}$$

The analytical solution of the infiltration equation along the *z*-axis can be expressed as:

$$h\_{\rm t}^{\*}(z,t) = \frac{2\left(1 - e^{\alpha h\_{\rm t}}\right)}{Lc} \mathbf{e}^{\alpha \cos \beta (L-z)/2} \sum\_{m=1}^{\infty} (-1)^{m} \left(\frac{\lambda\_{\rm m}}{\mu\_{\rm m}}\right) \sin(\lambda\_{\rm m} z) \mathbf{e}^{-\mu\_{\rm m} t} \tag{5.3}$$

**Fig. 5.1** Schematic diagram of a homogeneous unsaturated slope

108 5 Slope Stability Analysis Based on Analytical and Numerical Solutions

$$h\_s^\*(z) = \left(1 - \mathbf{e}^{\alpha h\_4}\right) \frac{1 - \mathbf{e}^{-a \cos \beta z}}{1 - \mathbf{e}^{-a \cos \beta L}}\tag{5.4}$$

This application example simulated transient infiltration in an unsaturated slope using the Chebyshev spectral method (CSM). The slope thickness (*L*) and angle (β) are assumed to be 10 m and 33°, respectively (Fig. 5.1). The initial condition was unchanged at *h*(*z*, *t* = 0) = −10 m. The upper and lower boundary conditions of the soil can be expressed as *h*(*z* = 0) = − 10 m and *h*(*z* = *L*) = 0.

For this example, this chapter selected sandy soil and silty loam. The experimental data (Brooks and Corey 1964; Lu and Likos 2004) and fitting curve of the two soils are shown in Fig. 5.2. The fitted parameters including *k*s, θs, θr, and α are listed in Table 5.1. The total simulation time was 2 h and *N* = 40.

The pressure head calculated using the CSM was introduced into the infinite slope stability analysis, and the factor of safety (*F*s) could be obtained (Iverson 2000; Liu et al. 2017):

$$F\_s = \frac{c'}{z\chi\_1 \cos \beta \sin \beta} + \frac{\tan \phi'}{\tan \beta} - \frac{h\chi\_\text{w} \tan \phi'}{z\chi\_1 \cos \beta \sin \beta} \tag{5.5}$$

where *c*' denotes the effective cohesion, φ' denotes the effective friction angle, γ<sup>t</sup> denotes the unit weight of soil, γw denotes the unit weight of water, and *z* denotes the soil thickness. It was assumed that the unit weight of soil was 19 kN/m3, and the effective cohesiveness and effective friction angle of sandy soil and silty loam were 0 kPa, 32° and 10 kPa, 20°, respectively (Lu and Likos 2004).

**Fig. 5.2** Fitting curves of the SWCC: **a** sandy soil and **b** silty loam


**Table 5.1** Parameters of unsaturated soils for analysis

### 5.2 Application of Unsaturated Soil Slope Stability Under Rainfall 109

The revised *F*s for analyzing slope stability can be obtained as (Liu et al. 2017):

$$F\_s = \frac{c'}{z\chi\_l \cos \beta \sin \beta} + \frac{\tan \phi'}{\tan \beta} - \frac{h\chi\_\mathbf{v} \left(\frac{\theta - \theta\_\mathbf{r}}{\theta\_\mathbf{v} - \theta\_\mathbf{r}}\right) \tan \phi'}{z\chi\_l \cos \beta \sin \beta} \tag{5.6}$$

The profiles of the pressure head of the two soils over time were compared with the analytical solutions. Figure 5.3 shows that the computed results were in good agreement with the analytical solutions. The numerical findings for sandy soil and silty loam indicate that the accuracy of the absolute error can reach the order of 10−<sup>5</sup> and 10−6, respectively (Fig. 5.4).

**Fig. 5.3** Results comparison with the exact solution for unsaturated soil: **a** sandy soil and **b** silty loam

**Fig. 5.4** Absolute error of the computed results with the analytical solutions for unsaturated soil: **a** sandy soil and **b** silty loam

Figure 5.5 demonstrates the computed results of *F*s for different unsaturated soils. Additionally, the range of *F*s for different soils is different, and *F*s for silty loam is significantly smaller than that for sandy soil. Table 5.2 shows variations in *F*s at different depths over the infiltration time for different soils. It can be seen that the factor of safety decreased with increasing depth and infiltration time. *F*s for sandy soil was greater than 1.0 in Table 5.2, and *F*s for silty loam was less than 1.0 at *z*  = 9.263 m. Consequently, the stability of the slope of sandy soil was better than that of the slope of silty loam. The comparison of the computed results reveals that rainfall-induced landslides are closely related to soil types.

**Fig. 5.5** Computed results of *F*s for unsaturated soil: **a** sandy soil and **b** silty loam



### *5.2.2 Application of P-CSIM in the Stability Analysis of Shallow Rainfall Slope*

The proposed P-CSIM is used to solve the Richards' equation for the slope, then the numerical solution of the pressure head is substituted into the infinite slope model, and the stability of the slope is analyzed. The factor of safety of the slopes is calculated as Eqs. (5.5)–(5.6).

### **5.2.2.1 Homogeneous Soil Slope**

The soil slope thickness is 150 cm, the slope angle is 33°, the initial condition for removing the boundary point is *h*(*z*, *t* = 0) = − 100 cm, and the boundary conditions are *h*(*z* = 0) = − 100 cm and *h*(*z* = *L*) = 0. The homogeneous unsaturated soil is assumed to be a sandy soil with a saturated permeability coefficient *k*s, residual volumetric water content θr, and saturated volumetric water content θs of 9 × 10−2 cm/h, 0.08, and 0.43, respectively (Liu et al. 2017). The desaturation coefficient α is 1 × 10−2. In addition, the unit weight, effective cohesiveness, and effective friction angle of the soil are 21.5 kN/m3, 4.6 kPa, and 30°, respectively (Lu and Likos 2004).

Figures 5.6 and 5.7 represent the computed profile of the factor of safety and the pressure head for the homogeneous soil slope. The suction (negative pressure head) in shallow parts of the slope decreases with the increase of rainfall time. In addition, it is found that *F*s decreases with the increase of rainfall time at the shallow slope zone. The unstable slope is mostly 130–150 cm deep, which may be strength deterioration and softening of the shallow layer caused by rainfall infiltration (Zhuang et al. 2018).

### **5.2.2.2 Two-Layer Soil Slopes**

For slopes of two-layer soils, the thickness of unsaturated Soil 1 (*L*1) and Soil 2 (*L*2) is assumed to be 130 cm and 20 cm, respectively (Fig. 5.8). The slope angle and boundary and initial conditions are the same as the slope for homogeneous soils. The physical parameter settings are listed in Table 5.3. It is assumed that Soil 1 is loess, and that the permeability coefficient of Soil 2 is greater than that of Soil 1.

The calculated results (Fig. 5.9) show that increase of the pressure head in layered unsaturated soils dominates slope stability. It can be found that the pressure head changes faster at the interface. Besides, because the permeability coefficient of the bottom layer is lower than that of the upper layer, the top layer (unsaturated Soil 2) may become unstable. Figure 5.10 demonstrates *F*s profile of the two-layer unsaturated soil slope. It shows that unstable slope ranges from the interface to the soil surface after *t* = 10 h. Figure 5.11 indicates that *F*s varies with duration at the interface (*z* = 130 cm) for unsaturated soil slopes. It is found that *F*s is strongly affected by duration, and that the upper layer of slope becomes unstable after *t* = 8.2 h.

**Fig. 5.6** Computed profile of pressure head for homogeneous unsaturated soil slopes

**Fig. 5.7** Computed profile of safety factor for homogeneous unsaturated soil slopes

**Fig. 5.8** Schematic diagram of two-layer soil slope model


**Table 5.3** Parameters used in two-layer soil slopes

Accordingly, the results illustrate that the pressure head caused by rainfall infiltration is closely related to the soil layer, and the interface of the soil layer plays a key role in the slope stability related to rainfall infiltration.

### *5.2.3 Application of Improved Picard Method into Unsaturated Slopes*

In this example, this chapter investigates 1D transient infiltration into the two-layer soil slopes. The van Genuchten model is adopted here. In Fig. 5.8, the thickness of soil layer 1 (*L*1) and soil layer 2 (*L*2) is 3 m and 2 m, respectively, and the slope angle (β) is 33°. The governing equation is shown in Eq. (5.1). The VGM model parameters of the two soils are listed in Table 5.4 (Zha et al. 2013). The bottom boundary condition is groundwater table, defined as *h*(*z* = 0, *t*) = 0. The top boundary condition is given by (Wu et al. 2020):

$$\left[K\frac{\partial h}{\partial z} + K\right]\_{z=L\_1+L\_2} = q\_l, \quad t < t\_l \tag{5.7}$$

**Fig. 5.9** Computed profile of pressure head for two-layer unsaturated soil slopes

**Fig. 5.10** Computed profile of safety factor for two-layer unsaturated soil slopes

**Fig. 5.11** Results of safety factor at the interface (*z* = 130 cm) for two-layer unsaturated soil slopes


**Table 5.4** Parameter of van Genuchten model

$$h\_{|z=L\_1+L\_2} = h\_\mathsf{I}, \quad \mathsf{t} \ge \mathsf{t}\_\mathsf{I} \tag{5.8}$$

where *t*I represents ponding time; *q*t denotes the rainfall flux at the soil surface when the rainfall time is less than the ponding time; and *h*t denotes constant pressure head after the ponding time. Here, *q*t is assumed to be *K*s2/4 and *h*t = 0. The initial condition is given by *h*(*z*, *t* = 0) = −*z* × 10 m.

The total simulation time and time step are taken to be 2 h and 0.01 h, respectively. Let the number of nodes be *N*1 = 60 and *N*2 = 40. Furthermore, NTG-PI is applied to solve Eq. (5.1), where μ = 1 is adopted. The unit weight of the soil, the effective cohesion, and friction angle are 19.5 kN/m3, 4.6 kPa, and 30°, respectively.

Figure 5.12 shows that the pressure head increases over duration, while the pressure head at the interface has a greater increase. Compared with NTG-PI, the pressure head obtained by traditional method PI has a greater increase at the interface as the rainfall time increases. This may underestimate the safety factor of the interface for two-layer soil slopes (Fig. 5.13). Moreover, the results illustrate that the soil layer structure will affect the distribution of the pressure head in the two-layer soil slopes during rainfall. When the hydraulic conductivity of the lower soil is lower than that of the upper one, it is easy to form ponding at the interface and induce landslides.

**Fig. 5.12** Computed profile of pressure head for two-layer soil slopes

**Fig. 5.13** Computed profile of *F*s for two-layer soil slopes


**Table 5.5** Parameters for different cases

### *5.2.4 Parameter Sensitivity Analysis of Slope Stability Under Rainfall*

The test uses the Gardner model, as described by Eqs. (1.13)–(1.15). The slope soil layer is assumed to be homogeneous soil, the mathematical model is depicted in Fig. 5.1, its thickness is assumed to be 2 m, and the saturated and residual water contents are set to 0.46 and 0.1. The distribution characteristics of saturated permeability coefficient, desaturation coefficient α, slope angle, and rainfall *q* on porewater pressure are examined, and their influence on soil slope stability is discussed. The parameter settings of different conditions are shown in Table 5.5. For the slope stability analysis, the soil unit weight, effective cohesion, and effective internal friction angle are 19.9 kN/m3, 10 kPa, and 26°, respectively. The upper boundary is set as the flow boundary, the lower boundary is the groundwater level, its pressure head is set to 0, and the initial condition for removing the boundary points is *h*(*z*, *t* = 0) = −*z* m.

### **5.2.4.1 Influence of Saturated Permeability Coefficient on Slope Stability**

Figure 5.14 depicts the effect of three different saturated permeability coefficients on the pressure head distribution in Case 1. With the increase of the saturated permeability coefficient of the soil, the distribution of the pressure head moves faster, and

**Fig. 5.14** Influence of different saturated permeability coefficients on pressure head distribution (Case 1)

with the increase of the rainfall time, the change of the pressure head is large. The numerical results illustrate that in the numerical simulation of rainfall infiltration for soil slope, the size of the saturated permeability coefficient can directly affect the speed of the pressure head distribution. Figure 5.15 indicates that the increase of the saturated permeability coefficient has a great influence on the stability of the shallow soil slope.

### **5.2.4.2 Influence of Desaturation Coefficient on Slope Stability**

Figure 5.16 depicts the effect of different desaturation coefficients on the pressure head distribution in Case 2. As the desaturation coefficient increases, the pressure head distribution lags behind. At the same time, the curvature of the pressure head distribution becomes obvious with the increase of the desaturation coefficient.

Figure 5.17 demonstrates that the slope safety factor decreases as the desaturation factor decreases.

### **5.2.4.3 Influence of Slope Angle on Slope Stability**

Figure 5.18 represents the effect of different slope angles on the pressure head distribution in Case 3. As the slope angle increases, the pressure head distribution shifts gradually to the right. At the same time, the curvature of the pressure head distribution becomes obvious with the increase of the slope angle. Figure 5.19 shows that the slope angle has an inherent key effect on the slope safety coefficient, and the slope safety factor decreases continuously with the increase of the slope angle.

### **5.2.4.4 Influence of Rainfall on Slope Stability**

Figure 5.20 depicts the effect of different rainfall on the pressure head distribution in Case 4. With the increase of rainfall, the migration speed of the pressure head distribution increases, which means that the soil matric suction dissipates faster. At the same time, the curvature of the pressure head distribution becomes obvious with the increase of rainfall. Figure 5.21 verifies that the slope safety factor decreases with increasing rainfall.

To sum up, in the analysis of rainfall slope stability, the slope angle has an inherent key role in the factor of safety, and the saturation permeability coefficient, desaturation coefficient, and rainfall all have a large influence on the distribution of pressure head. For the short-term rainfall period, when the desaturation coefficient is small and the saturated permeability coefficient is large, the larger the rainfall, the larger the rainwater infiltration depth, which will lead to the softening of the soil, the reduction

**Fig. 5.16** Influence of different desaturation coefficients on pressure head distribution (Case 2)

**Fig. 5.17** Relationship between the slope safety factor (*z* = 1 m) and the desaturation factor

**Fig. 5.18** Influence of different slope angles on pressure head distribution (Case 3)

**Fig. 5.19** Relationship between the slope safety factor (*z* = 0 m) and the slope angle

**Fig. 5.20** Influence of different rainfall on pressure head distribution (Case 4)

of the shear strength, and the shallow soil slope stability. As a result, factors such as desaturation coefficient, rainfall, permeability coefficient, and slope angle all affect the stability of unsaturated soil slopes.

### **5.3 Rainfall Landslide Case Study—Tung Chung Landslide**

According to the improved method, a comparative study is carried out with the monitoring data of the pore-water pressure of a slope in Tung Chung, Hong Kong, to verify the application effect of the proposed method. The long-term rainfall monitoring in this area provides systematic data for investigating rainfall-induced landslides. Many researchers have conducted extensive and in-depth studies on landslides in the Lantau area, including field experiments, statistical analysis, and remote sensing interpretation (Zhang et al. 2016). The mathematical model of this application case is shown in Fig. 5.1.

The measured data of rainfall and pressure head at monitoring point SP3 are shown in Fig. 5.22. The change of pressure head is basically consistent with the change of rainfall, which better reflects the change of infiltration boundary. The monitoring data can be divided into three periods according to the time of the rainfall peak. Among them, the calibration period is from 0 to 192 h (0:00 on June 8, 2001 to 0:00 on June 16, 2001), and the Bayesian stochastic inversion of soil hydraulic parameters can be performed according to this stage (Yang et al. 2018). The model parameters and permeability coefficients of the two soil–water characteristic curves are shown in Table 5.6. After calibrating the soil hydraulic parameters by random back analysis, the error between the model calculated and measured values is very small (Yang et al. 2018). Numerical simulations are carried out for two rainfall periods in the verification period and compared with the monitoring data. The first period is 5 h (from 4:00 on June 27 to 9:00 on June 27), and the average rainfall is 1.4 cm/h. The second period is 5 h (from 1:00 on July 7 to 6 on July 7), and the average rainfall is 2.2 cm/h. The third period is 3 h (23:00 on July 15 to 2:00 on July 16), and the average rainfall is 2.4 cm/h.

The water pressure meter SP3 is buried at a depth of 200 cm, and the groundwater table depth is 250 cm. The thickness of the soil layer (*L*) is 250 cm, and β is assumed to be 35°. The initial conditions and boundary conditions are expressed as follows:

$$\left[K\frac{\partial h}{\partial z} + K\cos\beta\right]\_{z=L} = q\tag{5.9}$$

$$h\_{|z=0} = 0\tag{5.10}$$

The space step is 2.5 cm, and the time step is 0.1 h. Figures 5.23 and 5.24 compare the results between the numerical solutions of the Richards' equation of the VGM and Gardner models and the measured values in the first period. It can be found that the Richards' equation described by the VGM model fits the measured values better. When *t* = 5 h, the relative error between the simulated and measured values is only 1.23%.

**Fig. 5.22** Measured data of rainfall and pressure head at monitoring point SP3


**Table 5.6** Soil–water characteristic curve and saturated permeability coefficient

**Fig. 5.23** Comparison between the numerical solution of Richards' equation for the VGM model and the measured value in the first period

**Fig. 5.24** Comparison between the numerical solution of Richards' equation for the Gardner model and the measured value in the first period

**Fig. 5.25** Comparison between the numerical solution of Richards' equation for the VGM model and the measured value in the second period

Figures 5.25 and 5.26 represent the comparison results between the numerical solutions of the Richards' equation of the VGM and Gardner models and the measured values in the second period. Compared with the Gardner model, the fitting effect of the Richards' equation described by the VGM model and the measured values is still better. The relative error is 9.22% when *t* = 5 h.

Figures 5.27 and 5.28 compare the numerical solutions of the Richards' equation and the measured values of the VGM and Gardner models in the third period. When *t* = 3 h, the relative error between the simulated value and the measured value computed by the two models is 3.54% and 3.67%, respectively. The numerical results of Richards' equation using the VGM model are closer to the measured values.

As shown in Table 5.7, in the three time periods, the overall root mean square error (RSE) and relative error (RE) of the VGM model and the measured data are both smaller than the values obtained by the Gardner model, which further indicates that the Richards' equation described by the VGM model has a better fitting effect with the measured data. Numerical results illustrate that the proposed method can well simulate the time-varying response of pressure head in the rainfall infiltration of unsaturated soil slopes and has a good application effect.

**Fig. 5.26** Comparison between the numerical solution of Richards' equation for the Gardner model and the measured value in the second period

**Fig. 5.27** Comparison between the numerical solution of Richards' equation for the VGM model and the measured value in the third period

**Fig. 5.28** Comparison between the numerical solution of Richards' equation for the Gardner model and the measured value in the third period


**Table 5.7** Error between the simulated and the measured values during the validation period

### **5.4 Conclusions**

This chapter studies the application of Chebyshev spectral method in slope stability analysis and the application of improved iterative methods P-CSIM and NTG-PI in shallow rainfall slope stability analysis. In order to verify the accuracy and practicability of the proposed improved method, a comparative study was carried out combining the monitoring data of slope in the Tung Chung area of Hong Kong and the Xiaoba landslide area. The conclusions are drawn as follows:

(1) The numerical results obtained by the CSM method are highly consistent with the transient analytical solution, which indicates that the proposed method is sufficiently accurate to handle transient infiltration problems associated with rainfall-induced landslides. The numerical solution obtained by the CSM method is introduced into the slope stability analysis to effectively evaluate the slope stability under rainfall conditions. The numerical solutions of the pressure head obtained by the improved methods P-CSIM and NTG-PI can effectively analyze the shallow landslides caused by rainfall. Combined with the rainfall data, a mathematical model was established for the monitoring points, the evolution process of the pressure head of the slope was solved by an improved iterative method, and the stability analysis of the unsaturated soil slope was carried out. The results indicate that the numerical solution is consistent with the measured value, and the relative error is small, showing a good application prospect. The numerical solution to coupled water flow and deformation in two-layer unsaturated porous media is obtained using a finite element method.

(2) A conceptual model of two-layer unsaturated soils is established to analyze the rainfall infiltration process under different conditions. A simplified analysis of an infinite slope is used to compute the factor of safety as a function of the depth of wetting front, and special attention is paid to the hydrological response at the interface between two layers of unsaturated porous media.

### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.