# Modeling of Transient Flow in Unsaturated Geomaterials for Rainfall-Induced Landslides Using a Novel Spacetime Collocation Method

Authors: Cheng-Yu Ku; Chih-Yu Liu; Yan Su; Jing-En Xiao
Journal: Geofluids
 (2018)

---

## Abstract

The modeling of transient flow in unsaturated soils for rainfall-induced landslides using a novel spacetime collocation method is presented. A numerical solution obtained in the spacetime coordinate system is approximated by superpositioning Trefftz basis functions satisfying the linearized Richards equation for collocation points on the spacetime domain boundary. The Gardner exponential model is adopted to derive the linearized Richards equation to describe the soil-water characteristic curve in unsaturated soils. To deal with the rainfall-induced landslides, the infinite slope stability analysis coupled with the proposed meshless method with the consideration of the fluctuation of time-dependent matric suction is developed. The proposed method is validated for several test problems. Application examples of transient modeling of flow for rainfall-induced landslides in homogenous unsaturated soils are also conducted. Numerical results demonstrate that the proposed method is highly accurate to deal with transient flow in unsaturated soils for rainfall-induced landslides. In addition, it is found that the numerical method using the Richards equation with the Gardner model may provide a promising solution for different soil textures.


---

## Body

## 1. Introduction

In recent decades, global climate change impacts have increased the frequency of severe meteorological phenomena such as heavy rainfall events [1, 2]. Increased precipitation events associated with extreme weather conditions affect the stability of slopes and can lead to a variety of geoenvironmental hazards [3, 4]. Consequently, shallow landslides may often occur within the vadose zone under unsaturated soil conditions [5–7]. The prediction of flow in transient conditions is important in engineering practice while considering the stability of unsaturated soil slopes [8–12] or poro-mechanical coupling influence for rainfall-induced shallow landslides in unsaturated soils [13, 14].In the past, theoretical models have been developed to deal with slope stability assuming a steady or quasi-steady water table associated with an infinite-slope stability analysis [15]. With regard to the influence of rainfall infiltration on slope stability in unsaturated soils, Iverson [16] assumed that nearly saturated soil is a simplified form of the linear diffusion equation of the Richards equation. However, the influence of matric suction of soil and the soil-water characteristic curve (SWCC) cannot be comprehensively considered during the analysis. The theory of geofluids when rainfall infiltrates vadose zones can be described using either the variably saturated flow equation or the generalized Richards equation [17, 18]. Governed by nonlinear physical relationships, the Richards equation is a highly nonlinear equation. The SWCC is a major factor that influences the nonlinear physical relationship of the hydraulic conductivity in unsaturated zones [19–22]. Since the Richards equation is highly nonlinear and an analytical solution can only be provided by some transformations [23], numerical solutions of the Richards equation are commonly used to model the flow process in unsaturated soils.Numerical approaches to the modeling of various unsaturated geofluid phenomena using the mesh-based methods, such as the finite difference method (FDM) or the finite element method (FEM), have been studied extensively [24, 25]. In spite of the success of using the conventional mesh-generated methods as effective numerical tools for the solution of transient flow problems in unsaturated soils, there remains growing interest in the development of novel advanced computational methods. Due to the flexibility in solving practical problems involving complex geometry and boundary conditions such as tidal fluctuation which is usually difficult for finding the analytical solution, meshless methods have attracted considerable attention recently [26, 27]. For such problems involving regions of irregular geometry, the use of numerical methods, particularly the boundary-type meshless method, to approximate numerical solutions is advantageous. In the past decades, several meshless methods have been proposed, such as the collocation Trefftz method (CTM) [28–31], the method of fundamental solutions (MFS) [32–34], smoothed particle hydrodynamics (SPH) [35], diffuse element method [36], and the boundary particle method [37]. Among these meshless methods, the CTM can be categorized into the boundary-type meshless method and provides the most accurate solutions with optimal numerical stability [38, 39]. Because the boundary-type meshless method is originally developed to deal with the boundary value problems in Euclidean space, the study of transient flow in unsaturated soils for rainfall-induced landslides using the proposed method has not been reported yet.In this paper, a novel spacetime collocation method for modeling of transient flow in unsaturated soils for rainfall-induced landslides is presented. For the modeling of transient flow in unsaturated soils, we proposed a pioneering work based on the proposed boundary-type meshless method and utilized the spacetime coordinate system instead of that in the original Euclidean space [10]. A numerical solution obtained in the spacetime coordinate system is approximated by superpositioning Trefftz basis functions satisfying the linearized Richards equation for collocation points on the spacetime domain boundary. To deal with the rainfall-induced landslides, the infinite slope stability analysis coupled with the proposed meshless method with the consideration of the fluctuation of time-dependent matric suction is developed. The validity of the proposed method is established for several problems by comparing the results with the exact solutions. Application examples of transient modeling of flow for rainfall-induced landslides in homogenous unsaturated soils are also conducted. The formulation of the proposed method is described as follows.


## 2. Modeling of Transient Flow in Unsaturated Soils

### 2.1. The Linearized Richards Equation for the Inclined Slope

The generalized Richards equation represents the movement of water in saturated and unsaturated soils, which can be expressed as(1)K∇2H=SsSw+Ch∂H∂t,where H represents the total groundwater head, K represents hydraulic conductivity, Ss represents the specific storage, Sw represents effective saturation, Ch represents specific capacity, ∇ represents gradient, and t represents time. If porous media are heterogeneous and anisotropic, (1) becomes
(2)∇K∇H=SsSw+Ch∂H∂t.This is also known as the variably saturated flow equation to describe the groundwater flow in saturated and unsaturated zones simultaneously [17, 18]. Equation (2) can be rewritten as the three-dimensional form
(3)∂∂xKxh∂H∂x+∂∂yKyh∂H∂y+∂∂zKzh∂H∂z=SsSw∂H∂t+Ch∂H∂t,where x points down the ground surface, y points to the tangent of the topographic contour passing through the origin, z is the vertical coordinate, normal to the xy plane, and Kxh, Kyh, and Kzh are the unsaturated hydraulic conductivity functions in lateral direction and the vertical direction, respectively. The original Richards equation is developed only for horizontal ground plane. To deal with the movement of water in an unsaturated inclined slope as shown in Figure 1, the Richards equation must be revised. The total head can be expressed as follows.
(4)H=h+E,where h is the pressure head and E is the elevation head. The elevation head can be presented as follows.
(5)E=xsinβ+zcosβ,where β is the angle of inclination. By substituting (4) and (5) into (3), we obtained
(6)∂∂xKxh∂h∂x−sinβ+∂∂yKyh∂h∂y+∂∂zKzh∂h∂z+cosβ=SsSw+Ch∂h∂t.Figure 1
The definition of the elevation head in an unsaturated inclined slope.Considering only the unsaturated flow, (6) can be rewritten as
(7)∂∂xKxh∂h∂x−sinβ+∂∂yKyh∂h∂y+∂∂zKzh∂h∂z+cosβ=Ch∂h∂t.The specific moisture capacity function can be defined by(8)Ch=∂θm∂h,where θm is the moisture content. Because Kxh, Kyh, and Kzh and Ch are functions of the pressure head, the Richards equation is nonlinear. To solve (7), three characteristic functions are required including the SWCC, the unsaturated hydraulic conductivity function, and the specific moisture capacity function.Considering an one-dimensional transient flow problem, the duration of the rainfall events is far shorter than the transmission time of pore water pressure in thex and y axes. This consideration can be realistic for the inclined slope, especially for the steep slope covered by a thin soil mantle. Besides, the numerical solution of the Richards equation is computationally expensive [40]. In this study, the Richards equation in one-dimensional vertical flow is considered as follows.
(9)∂∂zKzh∂h∂z+cosβ=∂θm∂t.To normalize the hydraulic conductivity of unsaturated soils with respect to its saturated hydraulic conductivity, we obtained(10)Kr=KzhKs,where Ks is the saturated hydraulic conductivity and Kr is the relative hydraulic conductivity which is a function of the pressure head. Equation (11) can be obtained by substituting (10) into (9).
(11)∂∂zKr∂h∂z+∂Kr∂zcosβ=1Ks∂θm∂t.The high nonlinearity of physical behavior of unsaturated soils is rooted in the SWCC. The SWCC represents the relationship between water content and matric suction in unsaturated soils, which is important for the seepage analysis in the vadose zone. Many mathematical equations have been proposed to describe the SWCC. The Gardner and van Genuchten models of SWCC are commonly found in analytical and numerical solutions for flow problems, respectively [19, 41]. The Gardner exponential model is mostly used to derive analytical solutions of unsaturated flow problems, whereas the van Genuchten model is used in numerical implementations to simulate the physical behavior of unsaturated soils. Gardner [41] proposed a simple one-parameter exponential model as
(12)Se=eαgh,where αg is the pore size distribution parameter, Se is the effective saturation defined by normalizing volumetric water content with its saturated and residual values as Se=θm−θr/θs−θr, θs and θr denote the saturated and residual water content, respectively. The volumetric water content can be expressed as follows.
(13)θm=θr+θs−θreαgh.For the porous media including gilat loam and sandy loam, the Gardner model exhibits good match with the van Genuchten model [42]. The relative hydraulic conductivity is modeled by the Gardner exponential model [43, 44] as
(14)Kr=eαgh.Adopting the Gardner exponential model, the linearized Richards equation for the one-dimensional transient Richards equation can be derived as follows [6, 10, 43, 45].
(15)∂2h¯∂z2+αgcosβ∂h¯∂z=c∂h¯∂t,where c=αgθs−θr/Ks, h¯ is the pressure head of the linearized Richards equation which can be defined as h¯=eαgh−eαghd, and hd is the pressure head when the soil is dry.


### 2.2. The Novel Spacetime Collocation Method

In this study, the proposed spacetime collocation method rooted from the CTM. The Trefftz method is first proposed by the German mathematician Erich Trefftz [28]. Later, the CTM [29–31, 46] is commonly used for the numerical solution of partial differential equations. Since the CTM is a boundary-type meshless method, it approximates the solutions using the T-complete basis functions where the solutions can be expressed as a linear combination of the T-complete basis functions automatically satisfying governing equations [38]. The CTM requires the evaluation of a coefficient for each term in the series. The determination of coefficients may be obtained by solving the unknown coefficients in the linear simultaneous equations which are accomplished by collocation imposing the boundary conditions at a number of points. The CTM used to deal with the Laplace equation in the two-dimensional cylindrical coordinate system. The numerical solution of the Laplace equation is approximated by superpositioning the T-complete basis functions satisfying the governing equation [39]. The T-complete basis functions are composed of a set of linearly independent functions which can be obtained by finding the general solutions to the Laplace equation in cylindrical coordinates. Since the CTM is originally developed to deal with the boundary value problems in the past, the application of the CTM for solving transient problems is hardly found [47].In this paper, we proposed a novel spacetime collocation method, named the spacetime collocation method (SCM) [10]. The numerical approximation was developed by using the general solution of the one-dimensional linearized Richards equation. Different from the conventional collocation method based on a set of unstructured points in space only, a pioneering work in this study uses the SCM by collocating the boundary points in the spacetime coordinate system such that the initial and boundary conditions can both be treated as boundary conditions on the spacetime domain boundary.We assume that time and the speed of light are absolute physical quantities that play the role of the independent variable such that the spacetime coordinate system is ann-dimensional space and one-dimensional of time [10]. For instance, the one-dimensional linearized Richards equation is one-dimensional in space and one-dimensional in time. In the past, the collocation points of the meshless method are only applied on the boundary in Euclidean space, as shown in Figure 2(a), and the conventional time-marching scheme is also adopted for time discretization. However, in this study, the spacetime domain is a one-dimensional space and one-dimensional of time which constructs a rectangular shape, as shown in Figure 2(b). Meanwhile, we transformed the one-dimensional initial value problem, into a two-dimensional inverse boundary value problem. It should be noted that the initial and boundary conditions are both applied on the spacetime boundary. In addition, it becomes an inverse boundary value problem because the final time boundary values may not be assigned.Figure 2
Illustration of the proposed collocation scheme.
(a)
Conventional collocation scheme for one-dimensional transient flow problem
(b)
Collocation points of one-dimensional transient flow problem in the spacetime domainConsidering the one-dimensional spacetime domain,Ωt, enclosed by a spacetime boundary, Γt, the one-dimensional linearized Richards equation can be expressed as
(16)∂2h¯∂z2+αgcosβ∂h¯∂z=c∂h¯∂tinΩt.The initial condition can be described as(17)h¯z,t=0=g¯,where g¯ denotes the distribution of the pressure head in the spacetime domain, Ωt, at time zero. To solve (16), the following boundary conditions must be given.
(18)h¯z,t=f¯onΓDt,(19)h¯nz,t=w¯onΓNt.where ΓDt denotes the spacetime boundary where the Dirichlet boundary condition is given, ΓNtdenotes the spacetime boundary where the Neumann boundary condition is given, f¯ denotes the Dirichlet boundary condition in the spacetime domain, and w¯ denotes the Neumann boundary condition in the spacetime domain.The transient pressure head of the linearized Richards equation [48] can be expressed as
(20)h¯=h¯t+h¯s,where h¯s and h¯t denote the steady-state and transient pressure head of the linearized Richards equation, respectively. The transient linearized Richards equation is determined by substituting (20) into (16), which gives
(21)∂2h¯t∂z2+αgcosβ∂h¯t∂z=c∂h¯t∂t.The process of the separation of variables may be used by having the transient solution as(22)h¯tz,t=Zzϖt.Substituting the above equation into (21) and dividing by Zzϖt gives
(23)ϖt∂2Zz∂z2+αgcosβϖt∂Zz∂z=cZz∂ϖt∂t,(24)1Zz∂2Zz∂z2+αgcosβ∂Zz∂z=cϖt∂ϖt∂t.Each term in the above equation must be a constant for a nonzero solution. Accordingly, we obtained the following equations.(25)1Zz∂2Zz∂z2+αgcosβ∂Zz∂z=η=cϖt∂ϖt∂t,(26)∂2Zz∂z2+αgcosβ∂Zz∂z−ηZz=0,(27)c∂ϖt∂t−ηϖt=0,where η is the arbitrary constant. This leads to the characteristic equations
(28)m12+αgcosβm1−η=0,(29)cm2−η=0,where m1=−αgcosβ±αgcosβ2+4η/2, and m2=η/c. The above equations are simple ordinary differential equations that have solutions,
(30)Zz=A1em1az+A2em1bz,(31)ϖt=em2t,where m1a=−αgcosβ+αgcosβ2+4η/2 and m1b=−αgcosβ−αgcosβ2+4η/2.The general solution forh¯tz,t becomes
(32)h¯tz,t=Zzϖt=A1em1az+A2em1bzem2t,where A1 and A2 are arbitrary constants to be determined. Rewriting the above equations in terms of sine and cosine series, we obtained the following equation:
(33)h¯tx≈∑k=1MpkWkx,where x=z,t, pk=p1kp2kp3kT, and Wkx=W1kW2kW3k. pk are unknown coefficients, Wkx are the Trefftz basis functions, T is the transpose of the matrix, k is the positive integer, M is the order of the basis function, W1k=e−αgcosβz/2−μktsinλkz, W2k=e−αgcosβz/2−μkt, W3k=e−αgcosβz/2−μktcosλkz, λk=πk/L, μk=αgcosβ2/4+λk2/c and L is the characteristic length. Equation (33) can be discretized at a number of collocated points on the spacetime boundary using the initial and boundary conditions, which can be expressed as follows.
(34)h¯tixi≈∑k=1MpkWkixi,where xi=zi,ti and Wkixi=W1kiW2kiW3ki.A system of simultaneous linear equations can then be obtained as(35)Aα=B.where A=W111W211W311⋯W3M1W112W212W312⋯W3M2⋮⋮⋮⋮⋮W11IW21IW31I⋯W3MI, α=p11p21p31⋮p3M, B=g¯1g¯2g¯3⋮g¯I, A is a matrix of the Trefftz basis functions with the size of I×3M, α with the size of 3M×1 is a vector of unknown coefficients, B with the size of I×1 is a vector of given values from the Dirichlet boundary condition at collocation points, I is the number of collocation points, 3M is the number of the terms related to the order of the basis function as depicted in (34), p11,p21,p31,…,p3M are unknown coefficients to be evaluated, and g¯1,g¯2,g¯3,…,g¯I are the values of the Dirichlet boundary condition. The boundary and source points are collocated in the spacetime domain. Solving the above simultaneous linear equations using the commercial program MATLAB backslash operator, we may obtain the set of unknown coefficients for the spacetime domain, Ωt. To obtain the pressure head results for the spacetime domain, the inner collocation points in Ωt must be placed. The pressure head, h¯tz,t, at inner collocation points can then be found by (35) using Α and α for the spacetime domain, Ωt. Once we obtained the transient pressure head of the linearized Richards equation, h¯tz,t, the transient solution can be obtained using the following two equations.
(36)h¯z,t=h¯tz,t+h¯sz,(37)h¯sz=1−eαghd1−e−αgcosβz1−e−αgcosβL,where h¯sz is the steady-state pressure head of the linearized Richards equation [45]. Finally, the transient solution can be obtained as follows.
(38)hz,t=1αglnh¯z,t+eαghd.


## 2.1. The Linearized Richards Equation for the Inclined Slope

The generalized Richards equation represents the movement of water in saturated and unsaturated soils, which can be expressed as(1)K∇2H=SsSw+Ch∂H∂t,where H represents the total groundwater head, K represents hydraulic conductivity, Ss represents the specific storage, Sw represents effective saturation, Ch represents specific capacity, ∇ represents gradient, and t represents time. If porous media are heterogeneous and anisotropic, (1) becomes
(2)∇K∇H=SsSw+Ch∂H∂t.This is also known as the variably saturated flow equation to describe the groundwater flow in saturated and unsaturated zones simultaneously [17, 18]. Equation (2) can be rewritten as the three-dimensional form
(3)∂∂xKxh∂H∂x+∂∂yKyh∂H∂y+∂∂zKzh∂H∂z=SsSw∂H∂t+Ch∂H∂t,where x points down the ground surface, y points to the tangent of the topographic contour passing through the origin, z is the vertical coordinate, normal to the xy plane, and Kxh, Kyh, and Kzh are the unsaturated hydraulic conductivity functions in lateral direction and the vertical direction, respectively. The original Richards equation is developed only for horizontal ground plane. To deal with the movement of water in an unsaturated inclined slope as shown in Figure 1, the Richards equation must be revised. The total head can be expressed as follows.
(4)H=h+E,where h is the pressure head and E is the elevation head. The elevation head can be presented as follows.
(5)E=xsinβ+zcosβ,where β is the angle of inclination. By substituting (4) and (5) into (3), we obtained
(6)∂∂xKxh∂h∂x−sinβ+∂∂yKyh∂h∂y+∂∂zKzh∂h∂z+cosβ=SsSw+Ch∂h∂t.Figure 1
The definition of the elevation head in an unsaturated inclined slope.Considering only the unsaturated flow, (6) can be rewritten as
(7)∂∂xKxh∂h∂x−sinβ+∂∂yKyh∂h∂y+∂∂zKzh∂h∂z+cosβ=Ch∂h∂t.The specific moisture capacity function can be defined by(8)Ch=∂θm∂h,where θm is the moisture content. Because Kxh, Kyh, and Kzh and Ch are functions of the pressure head, the Richards equation is nonlinear. To solve (7), three characteristic functions are required including the SWCC, the unsaturated hydraulic conductivity function, and the specific moisture capacity function.Considering an one-dimensional transient flow problem, the duration of the rainfall events is far shorter than the transmission time of pore water pressure in thex and y axes. This consideration can be realistic for the inclined slope, especially for the steep slope covered by a thin soil mantle. Besides, the numerical solution of the Richards equation is computationally expensive [40]. In this study, the Richards equation in one-dimensional vertical flow is considered as follows.
(9)∂∂zKzh∂h∂z+cosβ=∂θm∂t.To normalize the hydraulic conductivity of unsaturated soils with respect to its saturated hydraulic conductivity, we obtained(10)Kr=KzhKs,where Ks is the saturated hydraulic conductivity and Kr is the relative hydraulic conductivity which is a function of the pressure head. Equation (11) can be obtained by substituting (10) into (9).
(11)∂∂zKr∂h∂z+∂Kr∂zcosβ=1Ks∂θm∂t.The high nonlinearity of physical behavior of unsaturated soils is rooted in the SWCC. The SWCC represents the relationship between water content and matric suction in unsaturated soils, which is important for the seepage analysis in the vadose zone. Many mathematical equations have been proposed to describe the SWCC. The Gardner and van Genuchten models of SWCC are commonly found in analytical and numerical solutions for flow problems, respectively [19, 41]. The Gardner exponential model is mostly used to derive analytical solutions of unsaturated flow problems, whereas the van Genuchten model is used in numerical implementations to simulate the physical behavior of unsaturated soils. Gardner [41] proposed a simple one-parameter exponential model as
(12)Se=eαgh,where αg is the pore size distribution parameter, Se is the effective saturation defined by normalizing volumetric water content with its saturated and residual values as Se=θm−θr/θs−θr, θs and θr denote the saturated and residual water content, respectively. The volumetric water content can be expressed as follows.
(13)θm=θr+θs−θreαgh.For the porous media including gilat loam and sandy loam, the Gardner model exhibits good match with the van Genuchten model [42]. The relative hydraulic conductivity is modeled by the Gardner exponential model [43, 44] as
(14)Kr=eαgh.Adopting the Gardner exponential model, the linearized Richards equation for the one-dimensional transient Richards equation can be derived as follows [6, 10, 43, 45].
(15)∂2h¯∂z2+αgcosβ∂h¯∂z=c∂h¯∂t,where c=αgθs−θr/Ks, h¯ is the pressure head of the linearized Richards equation which can be defined as h¯=eαgh−eαghd, and hd is the pressure head when the soil is dry.


## 2.2. The Novel Spacetime Collocation Method

In this study, the proposed spacetime collocation method rooted from the CTM. The Trefftz method is first proposed by the German mathematician Erich Trefftz [28]. Later, the CTM [29–31, 46] is commonly used for the numerical solution of partial differential equations. Since the CTM is a boundary-type meshless method, it approximates the solutions using the T-complete basis functions where the solutions can be expressed as a linear combination of the T-complete basis functions automatically satisfying governing equations [38]. The CTM requires the evaluation of a coefficient for each term in the series. The determination of coefficients may be obtained by solving the unknown coefficients in the linear simultaneous equations which are accomplished by collocation imposing the boundary conditions at a number of points. The CTM used to deal with the Laplace equation in the two-dimensional cylindrical coordinate system. The numerical solution of the Laplace equation is approximated by superpositioning the T-complete basis functions satisfying the governing equation [39]. The T-complete basis functions are composed of a set of linearly independent functions which can be obtained by finding the general solutions to the Laplace equation in cylindrical coordinates. Since the CTM is originally developed to deal with the boundary value problems in the past, the application of the CTM for solving transient problems is hardly found [47].In this paper, we proposed a novel spacetime collocation method, named the spacetime collocation method (SCM) [10]. The numerical approximation was developed by using the general solution of the one-dimensional linearized Richards equation. Different from the conventional collocation method based on a set of unstructured points in space only, a pioneering work in this study uses the SCM by collocating the boundary points in the spacetime coordinate system such that the initial and boundary conditions can both be treated as boundary conditions on the spacetime domain boundary.We assume that time and the speed of light are absolute physical quantities that play the role of the independent variable such that the spacetime coordinate system is ann-dimensional space and one-dimensional of time [10]. For instance, the one-dimensional linearized Richards equation is one-dimensional in space and one-dimensional in time. In the past, the collocation points of the meshless method are only applied on the boundary in Euclidean space, as shown in Figure 2(a), and the conventional time-marching scheme is also adopted for time discretization. However, in this study, the spacetime domain is a one-dimensional space and one-dimensional of time which constructs a rectangular shape, as shown in Figure 2(b). Meanwhile, we transformed the one-dimensional initial value problem, into a two-dimensional inverse boundary value problem. It should be noted that the initial and boundary conditions are both applied on the spacetime boundary. In addition, it becomes an inverse boundary value problem because the final time boundary values may not be assigned.Figure 2
Illustration of the proposed collocation scheme.
(a)
Conventional collocation scheme for one-dimensional transient flow problem
(b)
Collocation points of one-dimensional transient flow problem in the spacetime domainConsidering the one-dimensional spacetime domain,Ωt, enclosed by a spacetime boundary, Γt, the one-dimensional linearized Richards equation can be expressed as
(16)∂2h¯∂z2+αgcosβ∂h¯∂z=c∂h¯∂tinΩt.The initial condition can be described as(17)h¯z,t=0=g¯,where g¯ denotes the distribution of the pressure head in the spacetime domain, Ωt, at time zero. To solve (16), the following boundary conditions must be given.
(18)h¯z,t=f¯onΓDt,(19)h¯nz,t=w¯onΓNt.where ΓDt denotes the spacetime boundary where the Dirichlet boundary condition is given, ΓNtdenotes the spacetime boundary where the Neumann boundary condition is given, f¯ denotes the Dirichlet boundary condition in the spacetime domain, and w¯ denotes the Neumann boundary condition in the spacetime domain.The transient pressure head of the linearized Richards equation [48] can be expressed as
(20)h¯=h¯t+h¯s,where h¯s and h¯t denote the steady-state and transient pressure head of the linearized Richards equation, respectively. The transient linearized Richards equation is determined by substituting (20) into (16), which gives
(21)∂2h¯t∂z2+αgcosβ∂h¯t∂z=c∂h¯t∂t.The process of the separation of variables may be used by having the transient solution as(22)h¯tz,t=Zzϖt.Substituting the above equation into (21) and dividing by Zzϖt gives
(23)ϖt∂2Zz∂z2+αgcosβϖt∂Zz∂z=cZz∂ϖt∂t,(24)1Zz∂2Zz∂z2+αgcosβ∂Zz∂z=cϖt∂ϖt∂t.Each term in the above equation must be a constant for a nonzero solution. Accordingly, we obtained the following equations.(25)1Zz∂2Zz∂z2+αgcosβ∂Zz∂z=η=cϖt∂ϖt∂t,(26)∂2Zz∂z2+αgcosβ∂Zz∂z−ηZz=0,(27)c∂ϖt∂t−ηϖt=0,where η is the arbitrary constant. This leads to the characteristic equations
(28)m12+αgcosβm1−η=0,(29)cm2−η=0,where m1=−αgcosβ±αgcosβ2+4η/2, and m2=η/c. The above equations are simple ordinary differential equations that have solutions,
(30)Zz=A1em1az+A2em1bz,(31)ϖt=em2t,where m1a=−αgcosβ+αgcosβ2+4η/2 and m1b=−αgcosβ−αgcosβ2+4η/2.The general solution forh¯tz,t becomes
(32)h¯tz,t=Zzϖt=A1em1az+A2em1bzem2t,where A1 and A2 are arbitrary constants to be determined. Rewriting the above equations in terms of sine and cosine series, we obtained the following equation:
(33)h¯tx≈∑k=1MpkWkx,where x=z,t, pk=p1kp2kp3kT, and Wkx=W1kW2kW3k. pk are unknown coefficients, Wkx are the Trefftz basis functions, T is the transpose of the matrix, k is the positive integer, M is the order of the basis function, W1k=e−αgcosβz/2−μktsinλkz, W2k=e−αgcosβz/2−μkt, W3k=e−αgcosβz/2−μktcosλkz, λk=πk/L, μk=αgcosβ2/4+λk2/c and L is the characteristic length. Equation (33) can be discretized at a number of collocated points on the spacetime boundary using the initial and boundary conditions, which can be expressed as follows.
(34)h¯tixi≈∑k=1MpkWkixi,where xi=zi,ti and Wkixi=W1kiW2kiW3ki.A system of simultaneous linear equations can then be obtained as(35)Aα=B.where A=W111W211W311⋯W3M1W112W212W312⋯W3M2⋮⋮⋮⋮⋮W11IW21IW31I⋯W3MI, α=p11p21p31⋮p3M, B=g¯1g¯2g¯3⋮g¯I, A is a matrix of the Trefftz basis functions with the size of I×3M, α with the size of 3M×1 is a vector of unknown coefficients, B with the size of I×1 is a vector of given values from the Dirichlet boundary condition at collocation points, I is the number of collocation points, 3M is the number of the terms related to the order of the basis function as depicted in (34), p11,p21,p31,…,p3M are unknown coefficients to be evaluated, and g¯1,g¯2,g¯3,…,g¯I are the values of the Dirichlet boundary condition. The boundary and source points are collocated in the spacetime domain. Solving the above simultaneous linear equations using the commercial program MATLAB backslash operator, we may obtain the set of unknown coefficients for the spacetime domain, Ωt. To obtain the pressure head results for the spacetime domain, the inner collocation points in Ωt must be placed. The pressure head, h¯tz,t, at inner collocation points can then be found by (35) using Α and α for the spacetime domain, Ωt. Once we obtained the transient pressure head of the linearized Richards equation, h¯tz,t, the transient solution can be obtained using the following two equations.
(36)h¯z,t=h¯tz,t+h¯sz,(37)h¯sz=1−eαghd1−e−αgcosβz1−e−αgcosβL,where h¯sz is the steady-state pressure head of the linearized Richards equation [45]. Finally, the transient solution can be obtained as follows.
(38)hz,t=1αglnh¯z,t+eαghd.


## 3. Model Verification and Comparison

### 3.1. Verification Example 1

The first verification example is an one-dimensional transient subsurface flow problem in homogenous unsaturated soil. The unsaturated soil is initially dry until water begins to infiltrate the soil. The ponding on the ground surface is then maintained holding the pressure head to zero. This is known as the one-dimensional Green-Ampt problem [49]. The height of the column of soil (L) is 10 (m). The saturated hydraulic conductivity (Ks), pore size distribution parameter (αg), saturated water content (θs), and residual water content (θr) for this example are 9×10−5 (m/h), 10−4, 0.50, and 0.11, respectively. The total simulation time (T) is ten hours (h). The governing equation can be expressed as follows.
(39)∂∂zKr∂h∂z+∂Kr∂z=1Ks∂θm∂t.The initial condition was the soil in dry condition wherehd=−105 (m). Thus,
(40)hz,t=0=hd.The boundary conditions of the top and bottom boundaries are the pressure head-type boundary condition, which can be expressed as follows.(41)hz=0,t=hd,(42)hz=L,t=0.The analytical transient solution of this example can be found as follows [45].
(43)h¯tz,t=21−eαghdLceαgL−z/2∑k=1M−1kλkμksinλkze−μtk.Theh¯z,t can be obtained using the following equations.
(44)h¯z,t=h¯tz,t+h¯sz,(45)h¯sz=1−eαghd1−e−αgz1−e−αgL.Finally, the transient solutionhz,t can be obtained as follows.
(46)hz,t=1αglnh¯z,t+eαghd.In this example, there is one-dimensional in space and one-dimensional in time. It is clear that the spacetime domain is a rectangular shape. We adopted 600 boundary collocation points and a source point. The number of boundary points in space domain and in time domain was considered to be 200 and 400, respectively. The Dirichlet boundary values were given on boundary collocation points which collocated on three sides of the domain except the right side which is the final time solution to be found. The order of the basis function depicted in (34) and the characteristic length for the transient analysis were 60 and 10, respectively. To view the results clearly, the profiles of the computed results on different times were selected to compare with the exact solution. To examine the accuracy of the proposed method, we collocated 2496 inner points uniformly placed inside the rectangular spacetime domain. This example was solved using both the proposed method and the finite difference method (FDM). For the FDM analysis, the spatial and temporal discretization for the problem must be considered separately. The central difference approximation was adopted for the spatial discretization. On the other hand, the implicit scheme is adopted for time discretization. The space interval and time interval were 0.02 (m) and 0.01 (h), respectively.The profiles of the computed results from the proposed method and the FDM on different times were then compared with the exact solution, as shown in Figure3. It is found that the computed results agreed well with the exact solution. Comparing the maximum absolute error, it is found that the maximum absolute errors are in the order of 10−6 and 10−1 for the proposed meshless method and the FDM, respectively, as shown in Figure 4. In this example, we proposed an innovative concept which is based on the spacetime coordinate system. Consequently, the transient flow problems can be transformed into the inverse boundary problem in which the boundary values of the final time are unknown. The advantages of the proposed method are that the highly accurate results can be obtained and the convergence problem can be avoided since the one-dimensional transient example is solved without using the traditional time-marching scheme.Figure 3
Comparison of results for solving one-dimensional Green-Ampt problem.Figure 4
Comparison of the absolute error of the computed results and the FDM.


### 3.2. Verification Example 2

The second example under investigation is the transient modeling of one-dimensional flow in a homogenous unsaturated inclined slope. A column of soil is initially dry until water begins to infiltrate the soil. A pool of water at ground surface is then maintained holding the pressure head to zero. The thickness of the soil is 10 (m) and the slope angle is 20 degrees, as shown in Figure5. The total simulation time is two hours. The governing equation can be expressed as (11).Figure 5
Schematic illustration of homogenous unsaturated inclined slope.The initial condition was the soil in dry condition maintained ashd=−10m. The boundary conditions at the top and bottom of the soil can be expressed as follows.
(47)hz=0,t=hd,(48)hz=L,t=0.The transient analytical solution can be expressed as follows [45].
(49)h¯tz,t=21−eαghdLceαgcosβL−z/2∑k=1M−1kλkμksinλkze−μtk.In this study, the modeling of one-dimensional transient flow in the homogenous unsaturated inclined slope is conducted. To validate the applicability of the proposed Gardner model, we selected two different soils which are sandy soil and silty loam for this example. The data sets from Brooks and Corey [50] and Lu and Likos [51] were adopted. These data sets cover a broad range of porous media, ranging from sandy soil to silty loam. The fitted parameters including saturated hydraulic conductivity, pore size distribution parameter, saturated water content, and residual water content are listed in Table 1 [51].Table 1
Unsaturated soil parameters used in the verification example 2.
Parameters for SWCC
Soil type
Sandy soil
Silty loam
Saturated hydraulic conductivity (m/h)
6.00×10−2
1.50×10−3
Pore size distribution parameter (1/m)
1.60×10−2
8.00×10−3
Saturated water content
0.50
0.46
Residual water content
0.11
0.14We adopted 375 boundary collocation points and a source point. The number of boundary points in space domain and in time domain was considered to be 125 and 250, respectively. The Dirichlet boundary values were given on boundary collocation points which collocated on three sides of the domain using the analytical solution for the problem. The order of the basis function depicted in (33) and the characteristic length for the transient analysis were 60 and 10, respectively. To obtain the computed results of the pressure head at different times, we collocated 2496 inner points uniformly placed inside the rectangular domain. The profiles of the numerical solution on different times were selected to compare with the analytical solution.The fitting curves of the SWCC for two different types of soil are illustrated in Figure6, and the computed results for sandy soil and silty loam in this investigation are shown in Figure 7. Figure 7 indicates that the computed results agree very well with the analytical solution. Results obtained demonstrate that the accuracy of the absolute error can be reached to the order of 10−6 and 10−4 for sandy soil and silty loam, as depicted in Figures 8(a) and 8(b), respectively.Figure 6
Fitting curves of the SWCC.
(a)
Sandy soil
(b)
Silty loamFigure 7
Result comparison with the exact solution for sandy soil and silty loam.
(a)
Sandy soil
(b)
Silty loamFigure 8
The absolute error of the computed results for sandy soil and silty loam.
(a)
Sandy soil
(b)
Silty loam


### 3.3. Verification Example 3

Since the proposed spacetime collocation method is advantageous in solving problems of irregular geometry, the third verification example is the two-dimensional transient subsurface flow problem in homogenous unsaturated soil enclosed by an irregular boundary. With a two-dimensional simply connected domain, the governing equation can be expressed as follows.(50)∂∂xKr∂h∂x+∂∂zKr∂h∂z+∂Kr∂zcosβ=1Ks∂θm∂t.The unsaturated soil parameters are the same as in verification example 1. The total simulation time is two hours. The initial condition is the unsaturated soil in dry condition maintained as pressure head to be –100 (m), which can be described as follows.(51)hx,z,t=0=−100.The two-dimensional irregular boundary under consideration can be defined as the following parametric equation.(52)Ω=x,z,tx=ρθcosθ,z=ρθsinθ.where ρθ=esinθ/2sin2θ2+esin2θcos2θ2/2,0≤θ≤2π.The boundary conditions are assumed to be the Dirichlet boundary condition. The Dirichlet boundary data is applied using the following exact solution [10].
(53)h¯tx,z,t=xe−αgcosβ2t/4ce−αgcosβz/2.The transient numerical solution can be obtained using the following equations [45].
(54)h¯x,z,t=h¯tx,z,t+h¯sx,z,(55)h¯sx,z=1−eαghdsinπxLeαgcosβ/2L−zsinhζzsinhζL.where ζ=αgcosβ2/4+π/L2. According to the previous study [10], the characteristic length is assumed to be 180. The CTM for the two-dimensional transient linearized Richards equation can also be composed of a set of linearly independent vectors adopting the method of the separation of variables. The solution can then be derived as follows.
(56)h¯tx,z,t≈∑k=1M∑j=1Msjkvjkx,z,t,where M is the order of the Trefftz basis function for approximating the solution. sjk are the unknown coefficients, which can be expressed as sjk=s1jks2jks3jks4jkT. s1jk, s2jk, s3jk, and s4jk are unknown coefficients, and vjkx,z,t are the Trefftz basis functions, which can be expressed as vjkx,z,t=v1jkv2jkv3jkv4jk, v1jk=e−αgcosβz/2−γjktsinλjxsinλkz, v2jk=e−αgcosβz/2−γjktsinλjxcosλkz, v3jk=e−αgcosβz/2−γjktcosλjxcosλkz, v4jk=e−αgcosβz/2−γjktcosλjxsinλkz, γjk=1/cζj2+λk2, ζj=αgcosβ2/4+λj2, and λj=πj/L.We utilized the spacetime coordinate system to model the transient two-dimensional flow in unsaturated soil. There is two-dimensional in space and one-dimensional in time in this two-dimensional transient analysis. Consequently, the spacetime domain is transformed into a three-dimensional irregular object domain, as shown in Figure9. Since the top side boundary values of the spacetime domain were unknowns to be determined, we transformed the two-dimensional initial value problem into the three-dimensional inverse boundary value problem. The initial and boundary conditions were applied on the bottom side and the circumferential irregular boundary of the spacetime domain, respectively.Figure 9
Spacetime boundary collocation points of the two-dimensional transient flow in unsaturated soils.The order of the basis function for this example was 10. We collocated 2162 boundary collocation points and a source point. The number of boundary points in the space domain and in the time domain was considered to be 162 and 2000, respectively. We adopted 1325 inner points uniformly placed inside the three-dimensional spacetime domain to obtain the computed results of the pressure head at different times. Figure10 shows the comparison of the computed result with the exact solution. It is found that the computed results agreed closely with the exact solution. Figure 11 depicts the absolute error of the two-dimensional computed results. Highly accurate numerical solutions in the order of 10−11 can be obtained for this two-dimensional transient analysis.Figure 10
Comparison of computed results with the exact solution for the modeling of two-dimensional transient flow in unsaturated soils.Figure 11
The absolute error of the computed results for the modeling of two-dimensional transient flow in unsaturated soils.


## 3.1. Verification Example 1

The first verification example is an one-dimensional transient subsurface flow problem in homogenous unsaturated soil. The unsaturated soil is initially dry until water begins to infiltrate the soil. The ponding on the ground surface is then maintained holding the pressure head to zero. This is known as the one-dimensional Green-Ampt problem [49]. The height of the column of soil (L) is 10 (m). The saturated hydraulic conductivity (Ks), pore size distribution parameter (αg), saturated water content (θs), and residual water content (θr) for this example are 9×10−5 (m/h), 10−4, 0.50, and 0.11, respectively. The total simulation time (T) is ten hours (h). The governing equation can be expressed as follows.
(39)∂∂zKr∂h∂z+∂Kr∂z=1Ks∂θm∂t.The initial condition was the soil in dry condition wherehd=−105 (m). Thus,
(40)hz,t=0=hd.The boundary conditions of the top and bottom boundaries are the pressure head-type boundary condition, which can be expressed as follows.(41)hz=0,t=hd,(42)hz=L,t=0.The analytical transient solution of this example can be found as follows [45].
(43)h¯tz,t=21−eαghdLceαgL−z/2∑k=1M−1kλkμksinλkze−μtk.Theh¯z,t can be obtained using the following equations.
(44)h¯z,t=h¯tz,t+h¯sz,(45)h¯sz=1−eαghd1−e−αgz1−e−αgL.Finally, the transient solutionhz,t can be obtained as follows.
(46)hz,t=1αglnh¯z,t+eαghd.In this example, there is one-dimensional in space and one-dimensional in time. It is clear that the spacetime domain is a rectangular shape. We adopted 600 boundary collocation points and a source point. The number of boundary points in space domain and in time domain was considered to be 200 and 400, respectively. The Dirichlet boundary values were given on boundary collocation points which collocated on three sides of the domain except the right side which is the final time solution to be found. The order of the basis function depicted in (34) and the characteristic length for the transient analysis were 60 and 10, respectively. To view the results clearly, the profiles of the computed results on different times were selected to compare with the exact solution. To examine the accuracy of the proposed method, we collocated 2496 inner points uniformly placed inside the rectangular spacetime domain. This example was solved using both the proposed method and the finite difference method (FDM). For the FDM analysis, the spatial and temporal discretization for the problem must be considered separately. The central difference approximation was adopted for the spatial discretization. On the other hand, the implicit scheme is adopted for time discretization. The space interval and time interval were 0.02 (m) and 0.01 (h), respectively.The profiles of the computed results from the proposed method and the FDM on different times were then compared with the exact solution, as shown in Figure3. It is found that the computed results agreed well with the exact solution. Comparing the maximum absolute error, it is found that the maximum absolute errors are in the order of 10−6 and 10−1 for the proposed meshless method and the FDM, respectively, as shown in Figure 4. In this example, we proposed an innovative concept which is based on the spacetime coordinate system. Consequently, the transient flow problems can be transformed into the inverse boundary problem in which the boundary values of the final time are unknown. The advantages of the proposed method are that the highly accurate results can be obtained and the convergence problem can be avoided since the one-dimensional transient example is solved without using the traditional time-marching scheme.Figure 3
Comparison of results for solving one-dimensional Green-Ampt problem.Figure 4
Comparison of the absolute error of the computed results and the FDM.


## 3.2. Verification Example 2

The second example under investigation is the transient modeling of one-dimensional flow in a homogenous unsaturated inclined slope. A column of soil is initially dry until water begins to infiltrate the soil. A pool of water at ground surface is then maintained holding the pressure head to zero. The thickness of the soil is 10 (m) and the slope angle is 20 degrees, as shown in Figure5. The total simulation time is two hours. The governing equation can be expressed as (11).Figure 5
Schematic illustration of homogenous unsaturated inclined slope.The initial condition was the soil in dry condition maintained ashd=−10m. The boundary conditions at the top and bottom of the soil can be expressed as follows.
(47)hz=0,t=hd,(48)hz=L,t=0.The transient analytical solution can be expressed as follows [45].
(49)h¯tz,t=21−eαghdLceαgcosβL−z/2∑k=1M−1kλkμksinλkze−μtk.In this study, the modeling of one-dimensional transient flow in the homogenous unsaturated inclined slope is conducted. To validate the applicability of the proposed Gardner model, we selected two different soils which are sandy soil and silty loam for this example. The data sets from Brooks and Corey [50] and Lu and Likos [51] were adopted. These data sets cover a broad range of porous media, ranging from sandy soil to silty loam. The fitted parameters including saturated hydraulic conductivity, pore size distribution parameter, saturated water content, and residual water content are listed in Table 1 [51].Table 1
Unsaturated soil parameters used in the verification example 2.
Parameters for SWCC
Soil type
Sandy soil
Silty loam
Saturated hydraulic conductivity (m/h)
6.00×10−2
1.50×10−3
Pore size distribution parameter (1/m)
1.60×10−2
8.00×10−3
Saturated water content
0.50
0.46
Residual water content
0.11
0.14We adopted 375 boundary collocation points and a source point. The number of boundary points in space domain and in time domain was considered to be 125 and 250, respectively. The Dirichlet boundary values were given on boundary collocation points which collocated on three sides of the domain using the analytical solution for the problem. The order of the basis function depicted in (33) and the characteristic length for the transient analysis were 60 and 10, respectively. To obtain the computed results of the pressure head at different times, we collocated 2496 inner points uniformly placed inside the rectangular domain. The profiles of the numerical solution on different times were selected to compare with the analytical solution.The fitting curves of the SWCC for two different types of soil are illustrated in Figure6, and the computed results for sandy soil and silty loam in this investigation are shown in Figure 7. Figure 7 indicates that the computed results agree very well with the analytical solution. Results obtained demonstrate that the accuracy of the absolute error can be reached to the order of 10−6 and 10−4 for sandy soil and silty loam, as depicted in Figures 8(a) and 8(b), respectively.Figure 6
Fitting curves of the SWCC.
(a)
Sandy soil
(b)
Silty loamFigure 7
Result comparison with the exact solution for sandy soil and silty loam.
(a)
Sandy soil
(b)
Silty loamFigure 8
The absolute error of the computed results for sandy soil and silty loam.
(a)
Sandy soil
(b)
Silty loam


## 3.3. Verification Example 3

Since the proposed spacetime collocation method is advantageous in solving problems of irregular geometry, the third verification example is the two-dimensional transient subsurface flow problem in homogenous unsaturated soil enclosed by an irregular boundary. With a two-dimensional simply connected domain, the governing equation can be expressed as follows.(50)∂∂xKr∂h∂x+∂∂zKr∂h∂z+∂Kr∂zcosβ=1Ks∂θm∂t.The unsaturated soil parameters are the same as in verification example 1. The total simulation time is two hours. The initial condition is the unsaturated soil in dry condition maintained as pressure head to be –100 (m), which can be described as follows.(51)hx,z,t=0=−100.The two-dimensional irregular boundary under consideration can be defined as the following parametric equation.(52)Ω=x,z,tx=ρθcosθ,z=ρθsinθ.where ρθ=esinθ/2sin2θ2+esin2θcos2θ2/2,0≤θ≤2π.The boundary conditions are assumed to be the Dirichlet boundary condition. The Dirichlet boundary data is applied using the following exact solution [10].
(53)h¯tx,z,t=xe−αgcosβ2t/4ce−αgcosβz/2.The transient numerical solution can be obtained using the following equations [45].
(54)h¯x,z,t=h¯tx,z,t+h¯sx,z,(55)h¯sx,z=1−eαghdsinπxLeαgcosβ/2L−zsinhζzsinhζL.where ζ=αgcosβ2/4+π/L2. According to the previous study [10], the characteristic length is assumed to be 180. The CTM for the two-dimensional transient linearized Richards equation can also be composed of a set of linearly independent vectors adopting the method of the separation of variables. The solution can then be derived as follows.
(56)h¯tx,z,t≈∑k=1M∑j=1Msjkvjkx,z,t,where M is the order of the Trefftz basis function for approximating the solution. sjk are the unknown coefficients, which can be expressed as sjk=s1jks2jks3jks4jkT. s1jk, s2jk, s3jk, and s4jk are unknown coefficients, and vjkx,z,t are the Trefftz basis functions, which can be expressed as vjkx,z,t=v1jkv2jkv3jkv4jk, v1jk=e−αgcosβz/2−γjktsinλjxsinλkz, v2jk=e−αgcosβz/2−γjktsinλjxcosλkz, v3jk=e−αgcosβz/2−γjktcosλjxcosλkz, v4jk=e−αgcosβz/2−γjktcosλjxsinλkz, γjk=1/cζj2+λk2, ζj=αgcosβ2/4+λj2, and λj=πj/L.We utilized the spacetime coordinate system to model the transient two-dimensional flow in unsaturated soil. There is two-dimensional in space and one-dimensional in time in this two-dimensional transient analysis. Consequently, the spacetime domain is transformed into a three-dimensional irregular object domain, as shown in Figure9. Since the top side boundary values of the spacetime domain were unknowns to be determined, we transformed the two-dimensional initial value problem into the three-dimensional inverse boundary value problem. The initial and boundary conditions were applied on the bottom side and the circumferential irregular boundary of the spacetime domain, respectively.Figure 9
Spacetime boundary collocation points of the two-dimensional transient flow in unsaturated soils.The order of the basis function for this example was 10. We collocated 2162 boundary collocation points and a source point. The number of boundary points in the space domain and in the time domain was considered to be 162 and 2000, respectively. We adopted 1325 inner points uniformly placed inside the three-dimensional spacetime domain to obtain the computed results of the pressure head at different times. Figure10 shows the comparison of the computed result with the exact solution. It is found that the computed results agreed closely with the exact solution. Figure 11 depicts the absolute error of the two-dimensional computed results. Highly accurate numerical solutions in the order of 10−11 can be obtained for this two-dimensional transient analysis.Figure 10
Comparison of computed results with the exact solution for the modeling of two-dimensional transient flow in unsaturated soils.Figure 11
The absolute error of the computed results for the modeling of two-dimensional transient flow in unsaturated soils.


## 4. Application Example

### 4.1. Hydrological Modeling of Geofluid in Unsaturated Soils

The application example under investigation is a one-dimensional, transient, unsaturated flow problem in a homogeneous inclined slope. In this example, the thickness of the soil is assumed to be 10 (m) and the slope angle is 33 degrees. The governing equation is depicted as (11). With regard to the influence of rainfall infiltration on slope stability in unsaturated zones, Iverson [16] proposed a one-dimensional linear diffusion equation which is a simplified form of the Richards equation. Accordingly, the influence of matric suction of soil and the SWCC cannot be comprehensively considered during the analysis using the one-dimensional linear diffusion equation. To demonstrate the difference of the result from this study and that from Iverson, we conducted a comparison example with the consideration of unsaturated flow in response to infiltration at ground surface. The total simulation time is two hours. Considering that the ponding on the ground surface is maintained holding the infiltration, the boundary conditions of the top and bottom boundaries are the pressure head-type boundary condition, which can be defined as (47) and (48). The soil was assumed to be glacial till where the following parameters are used in the example: Ks=1.8×10−3 (m/h), θs=0.32, and θr=0.03. The effective cohesiveness, effective friction angle, and unit weight of soil are 4.6 (kPa), 30 (degrees), and 21.5 (kN/m3), respectively. The unsaturated soil parameters used in this example are listed in Table 2 [51].Table 2
Unsaturated soil parameters used in the comparison example.
This study
van Genuchten model
Iverson’s model
Slope angle (degree)
33
Soli thickness (m)
10
Simulation time (h)
2
Saturated hydraulic conductivity (m/h)
1.8×10−3
Other parameters
αg=3.2×10−4
αv=2.8×10−3
D0=10~100×Ks
θs=0.32
nv=1.3
θr=0.03
θs=0.32
θr=0.03
Notation:αv is the fitting parameter related to the modal pore size, nv is the fitting parameter related to the pore size distribution, and D0 is the hydraulic diffusivity.To examine the applicability of the Gardner model, we compare the Gardner model with the van Genuchten model and measured data [20]. The van Genuchten parameters of unsaturated soil were obtained by fitting the following equation to the measured data.
(57)Se=1+αvhnv−mv,where αv is the fitting parameter related to the modal pore size, nv is the fitting parameter related to the pore size distribution, and mv=1−1/nv. The corresponding Gardner parameters were calculated using (12).From Figure12, it is found that the SWCC based on the Gardner model provides acceptable agreement with the van Genuchten model. Figure 13 reveals the computed results of the pressure head from the proposed method and those from van Genuchten and Iverson’s models. It is found that the computed results agree well with those from Iverson’s model, as depicted in Figures 13(a) and 13(b). However, the numerical solution from the simplified linear diffusion equation cannot obtain the appropriate pressure head when the value of the hydraulic diffusivity become smaller, as shown in Figures 13(c) and 13(d).Figure 12
Fitting results of the SWCC using van Genuchten and Gardner models.Figure 13
Comparison of computed pressure heads.
(a)
This study
(b)
van Genuchten model
(c)
Iverson’s model forD0=10×Ks
(d)
Iverson’s model forD0=100×Ks


### 4.2. Rainfall-Induced Slope Stability in Unsaturated Soils

To deal with the rainfall-induced landslides, the infinite slope stability analysis coupled with the proposed meshless method with the consideration of the fluctuation of time-dependent matric suction is developed. The failure mechanism of infinite slopes mainly involves colluvium, weathered rock formations, or shallow failure and planar slides in bedrock alternations located at thin depth below the ground level. The infinite slope theory is commonly used to evaluate the slope stability while the sliding soil thickness is far less than the slope height. The factor of safety (FS) for the infinite slope theory can be expressed as follows.
(58)FS=c'+zγtcos2β−hγwcos2βtanϕ'zγtcosβsinβ,where c' represents the effective cohesiveness, ϕ' represents the effective friction angle, z represents the soil thickness, γt represents the unit weight of soil, and γw represents the unit weight of water. Fredlund and Morgenstern adopted the perspectives of continuum mechanics and demonstrated that the shear strength equation can be applied to unsaturated soils [52]. Based on the SWCC, Vanapalli et al. proposed the shear strength equation suitable for unsaturated zones [53].
(59)τf=c'+σtanϕ'−hγwθm−θrθs−θrtanϕ',where σ represents the total stress and τf represents the shear strength. The revised shear strength equation was substituted into (58), and FS for analyzing slope stability can be obtained as follows.
(60)FS=c'zγtcosβsinβ+tanϕ'tanβ−hγwθm−θr/θs−θrtanϕ'zγtcosβsinβ.Figure14 demonstrates the computed results of FS. It is well known that Iverson’s model is based on the simplified Richards equation integrated with the infinite slope stability analysis for rainfall-induced shallow landside modeling [54]. Despite the fact that Iverson’s model may provide an effective tool for modeling landslides triggered by rainfall infiltration, the influence of matric suction and the high nonlinearity of physical behaviour of unsaturated soils cannot be considered comprehensively in Iverson’s model. The comparison of the computed results reveals that the numerical solution of the linear diffusion equation of the Richards equation can only be applied on the problems within a certain condition and should be used with caution when the value of the hydraulic diffusivity becomes smaller.Figure 14
Comparison of the computedFS.
(a)
This study
(b)
van Genuchten model
(c)
Iverson’s model forD0=10×Ks
(d)
Iverson model’s forD0=100×Ks


## 4.1. Hydrological Modeling of Geofluid in Unsaturated Soils

The application example under investigation is a one-dimensional, transient, unsaturated flow problem in a homogeneous inclined slope. In this example, the thickness of the soil is assumed to be 10 (m) and the slope angle is 33 degrees. The governing equation is depicted as (11). With regard to the influence of rainfall infiltration on slope stability in unsaturated zones, Iverson [16] proposed a one-dimensional linear diffusion equation which is a simplified form of the Richards equation. Accordingly, the influence of matric suction of soil and the SWCC cannot be comprehensively considered during the analysis using the one-dimensional linear diffusion equation. To demonstrate the difference of the result from this study and that from Iverson, we conducted a comparison example with the consideration of unsaturated flow in response to infiltration at ground surface. The total simulation time is two hours. Considering that the ponding on the ground surface is maintained holding the infiltration, the boundary conditions of the top and bottom boundaries are the pressure head-type boundary condition, which can be defined as (47) and (48). The soil was assumed to be glacial till where the following parameters are used in the example: Ks=1.8×10−3 (m/h), θs=0.32, and θr=0.03. The effective cohesiveness, effective friction angle, and unit weight of soil are 4.6 (kPa), 30 (degrees), and 21.5 (kN/m3), respectively. The unsaturated soil parameters used in this example are listed in Table 2 [51].Table 2
Unsaturated soil parameters used in the comparison example.
This study
van Genuchten model
Iverson’s model
Slope angle (degree)
33
Soli thickness (m)
10
Simulation time (h)
2
Saturated hydraulic conductivity (m/h)
1.8×10−3
Other parameters
αg=3.2×10−4
αv=2.8×10−3
D0=10~100×Ks
θs=0.32
nv=1.3
θr=0.03
θs=0.32
θr=0.03
Notation:αv is the fitting parameter related to the modal pore size, nv is the fitting parameter related to the pore size distribution, and D0 is the hydraulic diffusivity.To examine the applicability of the Gardner model, we compare the Gardner model with the van Genuchten model and measured data [20]. The van Genuchten parameters of unsaturated soil were obtained by fitting the following equation to the measured data.
(57)Se=1+αvhnv−mv,where αv is the fitting parameter related to the modal pore size, nv is the fitting parameter related to the pore size distribution, and mv=1−1/nv. The corresponding Gardner parameters were calculated using (12).From Figure12, it is found that the SWCC based on the Gardner model provides acceptable agreement with the van Genuchten model. Figure 13 reveals the computed results of the pressure head from the proposed method and those from van Genuchten and Iverson’s models. It is found that the computed results agree well with those from Iverson’s model, as depicted in Figures 13(a) and 13(b). However, the numerical solution from the simplified linear diffusion equation cannot obtain the appropriate pressure head when the value of the hydraulic diffusivity become smaller, as shown in Figures 13(c) and 13(d).Figure 12
Fitting results of the SWCC using van Genuchten and Gardner models.Figure 13
Comparison of computed pressure heads.
(a)
This study
(b)
van Genuchten model
(c)
Iverson’s model forD0=10×Ks
(d)
Iverson’s model forD0=100×Ks


## 4.2. Rainfall-Induced Slope Stability in Unsaturated Soils

To deal with the rainfall-induced landslides, the infinite slope stability analysis coupled with the proposed meshless method with the consideration of the fluctuation of time-dependent matric suction is developed. The failure mechanism of infinite slopes mainly involves colluvium, weathered rock formations, or shallow failure and planar slides in bedrock alternations located at thin depth below the ground level. The infinite slope theory is commonly used to evaluate the slope stability while the sliding soil thickness is far less than the slope height. The factor of safety (FS) for the infinite slope theory can be expressed as follows.
(58)FS=c'+zγtcos2β−hγwcos2βtanϕ'zγtcosβsinβ,where c' represents the effective cohesiveness, ϕ' represents the effective friction angle, z represents the soil thickness, γt represents the unit weight of soil, and γw represents the unit weight of water. Fredlund and Morgenstern adopted the perspectives of continuum mechanics and demonstrated that the shear strength equation can be applied to unsaturated soils [52]. Based on the SWCC, Vanapalli et al. proposed the shear strength equation suitable for unsaturated zones [53].
(59)τf=c'+σtanϕ'−hγwθm−θrθs−θrtanϕ',where σ represents the total stress and τf represents the shear strength. The revised shear strength equation was substituted into (58), and FS for analyzing slope stability can be obtained as follows.
(60)FS=c'zγtcosβsinβ+tanϕ'tanβ−hγwθm−θr/θs−θrtanϕ'zγtcosβsinβ.Figure14 demonstrates the computed results of FS. It is well known that Iverson’s model is based on the simplified Richards equation integrated with the infinite slope stability analysis for rainfall-induced shallow landside modeling [54]. Despite the fact that Iverson’s model may provide an effective tool for modeling landslides triggered by rainfall infiltration, the influence of matric suction and the high nonlinearity of physical behaviour of unsaturated soils cannot be considered comprehensively in Iverson’s model. The comparison of the computed results reveals that the numerical solution of the linear diffusion equation of the Richards equation can only be applied on the problems within a certain condition and should be used with caution when the value of the hydraulic diffusivity becomes smaller.Figure 14
Comparison of the computedFS.
(a)
This study
(b)
van Genuchten model
(c)
Iverson’s model forD0=10×Ks
(d)
Iverson model’s forD0=100×Ks


## 5. Conclusion

A novel spacetime collocation method for modeling of transient flow in unsaturated soils for rainfall-triggered landslides is developed. This pioneering work is based on the boundary-type meshless method and provides a promising solution for modeling the transient flow in unsaturated soils. The validity of the proposed method is established for a number of unsaturated flow problems. The fundamental concepts and the construct of the proposed method are clearly addressed in detail. Findings from this study are drawn as follows.(1)
The pioneering work in this study is the first successful attempt to solve the transient flow in unsaturated soils for an inclined slope using the novel spacetime collocation method. For the modeling of the transient flow in unsaturated soils, we proposed an innovative concept that one may utilize the spacetime coordinate system instead of that in the original Euclidean space. Consequently, both the initial and boundary conditions can be treated as boundary conditions on the spacetime domain boundary. The unsaturated flow problems can then be transformed into the inverse boundary value problem. Therefore, the one-dimensional problems can be solved without using the traditional time-marching scheme.
(2)
Since the SWCC is a major factor that influences the nonlinear physical relationship of unsaturated soils, we adopted the Gardner model to formulate the linearized Richards equation. For the numerical modeling of transient flow, we proposed a pioneering work using the proposed boundary-type meshless method in which the hydrological model coupled with the infinite slope stability analysis with the consideration of the SWCC can be used to model rainfall-induced landslides.
(3)
From the results of verification examples, it is found that the computed numerical results agree well with the analytical transient solution. Results from the numerical examples demonstrate that highly accurate numerical solutions with the error in the order of10−6 can be obtained. In addition, the proposed numerical scheme for the Richards equation to model one-dimensional flow in unsaturated inclined slopes is applicable for different soil textures.



---

*Source: 7892789-2018-08-08.xml*