# NUMERICAL MODELING OF UNSATURATED LAYERED SOIL FOR RAINFALL-INDUCED SHALLOW LANDSLIDES.

IntroductionDue to global climate change, severe weather phanomena such as extreme precipitation events are becoming much more frequent around the world (Pereira et al. 2010; Csete, Buzasi 2016). These affect the stability of slopes and have consequences on landslides. Since global warming is expected to increase the frequency and intensity of severe rainfall events, the number of people exposed to landslide risk may dramatically be increased (Gariano, Guzzetti 2016). Accordingly, it is of importance to understand and measure how climate variables and their variability affect landslide hazards. Landslides commonly observed in the roothill areas involve sliding surfaces which is typically called a shallow landslide with an indication of the high diffusion of shallow landslides phenomena in different environmental and geological contexts, all over the world. In shallow landslides, the sliding surface is located in the soil mantle or regolith typically to depths of a few centimeters to r few meters. Shallow landslides may often occur within the vadose zone under unsaturated soil conditions (Lu, Godt 2008; Yeh, Lee 2013). The assessment of rainfall-induced shallow landslides has long; been an interest of research topic in many branches of science and engineering such as soil science, hydrology, and geotechnical engineering (Montgomery Dietrich 1994; Baum et al. 20022; Frattini et al. 2004; Tsai, Yang 2006; Tsai, Chen 2010.

In the past, theoretical models have been developed to predict slope stability assuming a steady or quasi-steady water table associated with an infinite-slope stability analysis (Das 2010). With regard to the influence of rainfall infiltration on slope stability. in unsaturated zones, Iverson (2000) assumed that this occurred when soil was nearly saturated, further solving; the analytical solution of the Richards equation in the simplified form of a one-dimensional linear diffusion equation, and allowing the development of an analytical mod el of rainfall-induced shallow landslides. Baum et al. 12012) expanded Iverson's work and developed a transient rainfall infiltration and grid-based regional slope-stability analysis (TRIGRS) model. Because of its practicability, the TRIGRS model became popular for assessing shallow landslides (Crosta, Frattini 2003; Keim, Skaugset 2003; Lan et al. 2005; D'Odorico et al. 2005; Liu et al. 2008; Baum et al. 2010; Bordoni et al. 2015; Schiliro et al. 2015; Alvioli, Baum 2016).

In these numerical models, the hillslope soil is typically assumed to be homogeneous. Besides, Iverson's method assumed that nearly saturated soil was a simplified form of the linear diffusion equation of the Richards equation, in which the influence of matric suction of soil and the soil water characteristic curve (SWCC) cannot be comprehensively considered during analysis. A complete theory of groundwater flow when rainfall infiltrates unsaturated zones can be described using either the variably saturated flow equation or the generalized Richards equation. The Richards equation is highly nonlinear and cannot directly provide an analytical solution. The soil water content is a major factor that influences the nonlinear physical relationship of the hydraulic conductivity in unsaturated zones (Van Genuchten 1980; Fredlund et al. 1994).

Since the appearance of layered soil is much more common than homogeneous soil in nature, the hydrological process in heterogeneous porous media has drawn much attention and been studied (Fok 1970; Aylor, Parlange 1973; Hachum, Alfaro 1980; Samani et al. 1989), using analytical or numerical methods (Hanks Bowers 1962; Whisler, Klute 1966; Romano et al. 1988; Moldrup et al. 1989; Srivastava, Yeh 1991; Corradini et al. 2000; Ku, Tsai 2013). Because the numerical modeling of rainfall-induced shallow landslides in unsaturated layered soil using the variably saturated flow equation has hardly been reported, the study proposed a pioneer work on numerical modeling of rainfall-induced shallow landslides in unsaturated layered soil.

To model the shallow landslides, the infinite slope stability analysis coupled with the hydrological model with the consideration of the fluctuation of time-dependent pore water pressure and the SWCC proposed by Gardner was developed. A linearization process for the nonlinear Richards equation to deal with groundwater flow in unsaturated layered soil is derived using the Gardner model. To solve one-dimensional flow in the unsaturated zone of layered soil profiles, flux conservation and the continuity of pressure potential at the interface between two consecutive layers are considered in the numerical discretization of the finite difference method. In addition, model verification and comparison were performed, and a numerical model of rainfall-induced shallow landslides in unsaturated zones was subsequently developed. The formulation of the proposed method is described as follows.

1. Methods

1.1. Formulation of the variably saturated flow equation

Unsaturated soil is composed through a three-phase mixing process that involves a gas phase (air), a liquid phase (pore water) and a solid phase (soil particles). Soil layers above the water table are termed as the unsaturated zone. In geotechnical engineering, pore water pressure fluctuation reduces soil suction which may cause soil failure. Generally, continuous rainfall reduces soil suction, affecting the stability of slope. The complete theory of groundwater flow and rainfall infiltration of the unsaturated zone can be expressed using the Richards equation:

K [[nabla].sup.2]H = [[S.sub.s][S.sub.w] + C (h)][partial derivative]H/[partial derivative]t, (1)

where H represents the total groundwater head, K represents hydraulic conductivity, [S.sub.s] represents the specific storage, [S.sub.w] represents saturation, C(h) represents specific capacity, and t represents time. If porous media are heterogeneous and anisotropic, Eq. (1) becomes:

[nabla]K[nabla]H = [[S.sub.s][S.sub.w] + C (h)][partial derivative]H/[partial derivative]t. (2)

This is known as the generalized Richards equation or variably saturated flow equation, which can be used to describe the groundwater flow in saturated and unsaturated zones simultaneously (Tothova et al. 2007). Eq. (2) can be rewritten as follows for a complete three-dimensional problem:

[mathematical expression not reproducible],

where [K.sub.x] (h), [K.sub.y] (h), and [K.sub.z] (h) represent hydraulic conductivity in the x, y and z axis, which are all functions of the pore water pressure head h when unsaturated. To convert Eq. (3) into an equation that considers infiltration, we assume that the x and y planes represent slope vertical sections, the z axis is perpendicular to the xy plan, and the total water head H can be expressed using the pressure head h and elevation head E as follows:

H = h + E. (4)

When the slope infiltration problem is considered, the elevation head E varies with the slope angle a, as shown in Figure 1. The elevation head can be presented as follows:

E = x sin ([alpha]) + z cos ([alpha]). (5)

The governing equation can be obtained by substituting Eqs (4) and (5) into Eq. (3):

[mathematical expression not reproducible]. (6)

Thus, the majority of the variably saturated flow problem was simplified using Eq. (6). Considering a one-dimensional problem for example, rainfall duration was far shorter than the transmission time of pore water pressure in the x and y axis. This consideration can be realistic for slopes, especially steep slopes usually covered by a thin soil mantle. Thus, only vertical water flow variations were considered, which can be presented as a one-dimensional vertical infiltration equation:

[mathematical expression not reproducible]. (7)

When soil is saturated, [S.sub.w] = 1 and C (h) = 1. Therefore, the saturated flow equation is a commonly used governing equation for groundwater in the saturated zone, where:

[S.sub.s] [partial derivative]h/[partial derivative]t = [partial derivative]/[partial derivative]z] ([K.sub.z](h)([partial derivative]h/[partial derivative]z + cos [alpha])). (8)

When soil is not saturated, the groundwater flow equation in unsaturated zone can be expressed as:

C(h)[partial derivative]h/[partial derivative]t = [partial derivative]/[partial derivative]z] ([K.sub.z](h)([partial derivative]h/[partial derivative]z + cos [alpha])). (9)

When the soil layer is not saturated, Eq. (9) can be used to calculate fluctuations in pore water pressure. When the soil layer is gradually saturated by rainfall, Eq. (8) is used for analysis because Eq. (7) primarily describes the pore water pressure changes in saturated and unsaturated zones. This study adopted Eq. (7) as the core governing equation of rainfall-induced pore water pressure fluctuations in unsaturated zones. In 1958, Gardner (Gardner 1958) proposed a simple-one parameter model as follows:

[mathematical expression not reproducible], (10)

where [[alpha].sub.g] is the soil pore-size distribution parameter which is related to the pore size distribution of soil. [S.sub.e] is the effective saturation defined by normalizing volumetric water content with its saturated and residual values as:

[S.sub.e] = ([theta] - [[theta].sub.r]]/[[theta].sub.s] - [[theta].sub.r]), (11)

where [theta] represents the volumetric water content function, [[theta].sub.r] represents the residual water content and [[theta].sub.s] represents the saturated water content. Substituting Eq. (11) into Eq. (10), we have

[mathematical expression not reproducible]. (12)

It is common to normalize the hydraulic conductivity of unsaturated soil with respect to their maximum value. The normalized value, referred to as the relative hydraulic conductivity, [K.sub.r] can be expressed as:

[K.sub.r](h) = [K.sub.z](h)/[K.sub.s], (13)

where [K.sub.s] is the saturated hydraulic conductivity. [K.sub.r] (h) is the relative hydraulic conductivity and it is a function of the pressure head which makes the Richards equation highly nonlinear. Rather than using the van Genuchten expression for [S.sub.e], a simpler version is used as follows:

[S.sub.e] = [K.sub.r]. (14)

Therefore, [mathematical expression not reproducible]. Using the Gardner model, we can derive a linearized Richards equation. First of all, we define a new parameter [h.sup.*] which can be expressed as

[mathematical expression not reproducible], (15)

where [h.sup.*] is the pressure head of the linearized Richards equation and [chi] is a constant which is a key transformation allowing the solution of the linearized Richards equation using the transform methods (Duffy 1994; Tracy 2006, 2011; Liu et al. 2015). The parameter x has no physical meaning and can be expressed as [mathematical expression not reproducible]. [h.sub.d] is the pressure head when the soil is dry. Taking the derivative of Eq. (15) with respect to z, we obtain

[mathematical expression not reproducible]. (16)

Rearranging Eq. (16), we have

[mathematical expression not reproducible]. (17)

Again, substituting Eq. (15) into Eq. (17), we have

[mathematical expression not reproducible]. (18)

Eq. (18) can be written as

[K.sub.r] [partial derivative]h/[partial derivative]z = 1/[[alpha].sub.g] [partial derivative][h.sup.*]/[partial derivative]z. (19)

Taking the derivative of [mathematical expression not reproducible] with respect to z, we have

[mathematical expression not reproducible]. (20)

Substituting Eq. (17) into Eq. (20), we have

[partial derivative][K.sub.r]/[partial derivative]z = [partial derivative][h.sup.*]/[partial derivative]z. (21)

Recall Eq. (12), we have

[partial derivative][theta]/[partial derivative]t = ([[theta].sub.s] - [[theta].sub.r]) [partial derivative][S.sub.e]/[partial derivative]t. (22)

Because [mathematical expression not reproducible], we can rewrite the above equation as

[partial derivative][theta]/[partial derivative]t = ([[theta].sub.s] - [theta]r) [partial derivative][h.sup.*]/[partial derivative]t. (23)

Finally, we substitute Eqs. (17), (12), and (23) to Eq. (19). The linearized Richards equation can be obtained as

[mathematical expression not reproducible]. (24)

For simplicity, the equation can be expressed as

[K.sub.[alpha][theta]] [[partial derivative].sup.2][h.sup.*]/[partial derivative][z.sup.2] + [K.sub.[theta]] [partial derivative][h.sup.*]/[partial derivative]z = [partial derivative][h.sup.*]/[partial derivative]t, (25)

where

[K.sub.[theta]] = [K.sub.s]/([[theta].sub.s] - [[theta].sub.r]) and [K.sub.[alpha][theta]] = [K.sub.[theta]]/[[alpha].sub.g]. (26)

1.2. The finite difference method

The finite difference method (FDM) is adopted for the numerical discretization. The reason for choosing the FDM is that the FDM can be easily used to integrate the formulation for the discontinuities (or soil interface) which characterize the layers for unsaturated layered soil. In addition, several studies of numerical modeling of shallow landslides have been carried out using the FDM such as Yoshimatsu, Abe (2006), Wang et al. (2006), Singh et al. (2010), Sarkar et al. (2012). The main advantage of the FDM is that it is probably the most straight forward numerical method to formulate complicated partial different equations such as the Richards equation. In this study, we need to derive the sophisticated mathematical formulation for the unsaturated layered soil. Accordingly, we adopted the FDM as the numerical method. The linearized Richards equation based on the Gardner model is depicted in Eq. (25). To apply the FDM, we have to consider the discretization in spatial and temporal domains. For the spatial discretization, we divide the domain into sections, each of length Dz along the z axis and approximated the first and second derivatives in the linearized Richards equation for each grid point by the central difference formulas. For the time discretization, the implicit scheme is adopted. Accordingly, we obtain the finite difference equation as:

[mathematical expression not reproducible] (27)

where i is the nodal point and [DELTA]z is the interval. [DELTA]t is the time interval. [K.sub.s] [[alpha].sub.g], [[theta].sub.s] and [[theta].sub.r] are the same definition as described in the previous section. Considering the steady-state situation, we obtain the finite difference equation of the linearized Richards equation as:

[mathematical expression not reproducible], (28)

Eq. (28) can be used to deal with homogenous soil. For unsaturated layered soil, it is necessary to account for the discontinuities of the parameters which characterize the layers. This gives rise to the so-called interface problem for the linearized Richards equation. For the interface between layers, the conservation of water flux and the Darcy's law have to be considered, as shown in Figure 2. The finite difference equation of the first term, [K.sub.[alpha][theta]][[partial derivative].sup.2][h.sup.*]/[partial derivative][z.sup.2] Eq. (25), can be derived as follows:

[mathematical expression not reproducible], (29)

where [K.sub.[alpha][theta]A] and [K.sub.[alpha][theta]B] represent the parameter, [K.sub.[alpha][theta]], in layer A and layer B respectively. The finite difference equation of the second term, [K.sub.[theta]][partial derivative][h.sup.*]/[partial derivative]z in Eq. (25), can be derived as follows:

[mathematical expression not reproducible], (30)

where [K.sub.[theta]A] and [K.sub.[theta]B]- represent the parameter, [K.sub.[theta]], in layer A and layer B respectively. Substituting Eqs. (29) and (30) into Eq. (27) we obtain:

[mathematical expression not reproducible], (31)

where [mathematical expression not reproducible].

In Eqs. (28) and (31), [K.sub.[alpha][theta]] and [K.sub.[theta]] at the interface have to be determined. In a stratified soil deposit where the hydraulic conductivity for flow in a given direction changes from layer to layer, an equivalent hydraulic conductivity can be obtained by the following equation:

[K.sub.(eqv)] = H/([H.sub.1]/[K.sub.1]) + ([H.sub.2]/[K.sub.2]) + ... ([H.sub.n]/[K.sub.n]), (32)

where [K.sub.(eqv)] is the equivalent hydraulic conductivity, H is the height of the soil, [H.sub.i] is the thickness of each soil layer, and [K.sub.i] is the hydraulic conductivity of each soil layer. i is the number of the soil layer. Because we considered only two consecutive layers at the interface in the finite difference discretization. Using Eq. (32), [K.sub.[alpha][theta]] and [K.sub.[theta]] at the interface can be obtained as:

[mathematical expression not reproducible]. (33)

Accordingly, the finite difference equations of the linearized Richards equation for steady-state and transient conditions can be expressed as:

[mathematical expression not reproducible]; (34)

[mathematical expression not reproducible]. (35)

Furthermore, Eqs (34) and (35) can be expressed as the matrix form and can be written as

Ax = b, (36)

where A is a n x n square matrix, x is a n x 1 vector in which x = [[[h.sup.*.sub.1] ... [h.sup.*.sub.n]].sup.T], and b is a n x 1 vector. n is the number of unknowns to be solved. Eq. (36) can then be simply solved by x = [A.sup.-1]b .

1.3. Slope stability theory of unsaturated soil

The failure mechanism of infinite slopes primarily involves colluvium, weathered rock formations, or shallow failure and planar slides in bedrock alternations located at thin depth below the ground level. The sliding soil thickness is far less than the slope height. In most studies, infinite slopes are used to analyze shallow landslides. The factor of safety (FS) for the infinite slope theory can be expressed as follows:

FS = c' + (z[[gamma].sub.t] [cos.sup.2] [alpha] - h[[gamma].sub.w] [cos.sup.2] a[alpha]) tan [phi]' z[[gamma].sub.t] [cos.sup.2] [alpha] sin [alpha], (37)

where c' represents the effective cohesiveness, [phi]' represents the effective friction angle, z represents the soil thickness, [[gamma].sub.t] represents the unit weight of soil and [[gamma].sub.w] represents the unit weight of water. Fredlund adopted the perspectives of continuum mechanics and demonstrated that the shear strength equation can be applied to unsaturated soil (Fredlund, Morgenstern 1977). The reason of choosing Fredlund's approach for modeling the FS in unsaturated condition is that its description of the stress state of the unsaturated soil examine the context of multiphase continuum mechanics. Unlike saturated condition, which are two-phase mixing process involved a solid phase and a gas phase or a liquid phase, unsaturated soil is composed through a three-phase mixing process that comprises of a gas phase, a liquid phase and a solid phase. Therefore, the corresponding pressure of a gas and a liquid phase in unsaturated condition have a direct impact on physical behavior of the shear strength (Lu, Likos 2004). Vanapalli et al. proposed the shear strength equation based on the SWCC, which is suitable for unsaturated zones (Vanapalli et al. 1996):

[[tau].sub.f] = c' + [sigma] tan[phi]' - h[[gamma].sub.w] ([theta]-[[theta].sub.r]/[[theta].sub.s] - [[theta].sub.r]) tan [phi]', (38)

where [sigma] represents the total stress, [[tau].sub.f] represents the shear strength. The shear strength of unsaturated soil was developed to describe the physical consideration to the air-water interphase and these stress state variables (Fredlund et al. 1978; Vanapalli et al. 1996). The revised shear strength equation was substituted into the original infinite slope equation to obtain the following FS equation for analyzing slope stability in shallow unsaturated zones:

[mathematical expression not reproducible]. (39)

2. Verification examples

2.1. Verification example 1

The first example is a one-dimensional transient unsaturated groundwater flow problem in a homogenous soil. A column of soil is initially dry until water begins to infiltrate the soil. A pool of water at the ground surface is then maintained holding the pressure head to zero. The thickness of the soil L is 10 m. The initial condition is described as follows:

h (z, t = 0) = -[10.sup.5]. (40)

The boundary conditions of the top and bottom boundaries are as follows:

h (z = 0, t) = -[10.sup.5], (41)

h(z = L,t) = 0. (42)

Using Eq. (15), the suction head can be converted to [h.sup.*] as follows:

[h.sup.*] (z = 0, t) = 0, (43)

[h.sup.*] (z = L, t) = 1. (44)

In this example, the soil is assumed to be silt which the [[alpha].sub.g] of the Gardner model is[ 10.sup.-4]. The saturated hydraulic conductivity ([K.sub.s]), saturated water content ([[theta].sub.s]), and residual water content ([[theta].sub.r]) for this example are 9 x [10.sup.-5] m/s, 0.50, and 0.11, respectively, as shown in Table 1. For one-dimensional transient unsaturated groundwater flow in a homogenous soil, Tracy proposed the analytical transient solution as follows (Tracy 2011):

[mathematical expression not reproducible], (45)

where [[lambda].sub.k] = k[pi]/L, c = [[alpha].sub.g] ([[theta].sub.s] - [[theta].sub.r])/[K.sub.s], [mu]k = ([[alpha].sup.2.sub.g]/4 + [[lambda].sup.2.sub.k])/c, t is time. L is the thickness of the soil. Figure 3 demonstrates that the numerical results agreed well with the analytical solution for the one-dimensional transient unsaturated groundwater flow in a homogenous soil. It is found that the mean absolute error (MAE) was 1.18 x [10.sup.-4] m. The above verification example demonstrates that the proposed numerical scheme for the Richards equation to solve one-dimensional flow processes in the unsaturated zone of soil is validated.

2.2. Verification example 2

The second example under investigation is also a one-dimensional, transient, unsaturated groundwater flow problem in homogeneous soil. A soil column is initially dry until water begins to infiltrate the soil. A pool of water at the ground surface is maintained, holding the infiltration rate to zero. This is known as one-dimensional Green-Ampt problem (Green, Ampt 1911). The thickness of the soil L is 10 m. The initial condition is described as follows:

h (z, t = 0) = -100. (46)

The boundary conditions of the top and bottom boundaries are as follows.

h (z = 0, t) = -100, (47)

h (z = L,t) = 0. (48)

The soil was assumed to be sand, which has saturated hydraulic conductivity ([K.sub.s]) of 3 x [10.sup.-3] m/s. The saturated water content ([[theta].sub.s]) and residual water content ([[theta].sub.r]) are 0.50 and 0.11, respectively. In the Gardner model, there is only one fitting parameter, [[alpha].sub.g]. The SWCC obtained from the laboratory is shown in Figure 4. Using the Gardner model, [[alpha].sub.g] can be obtained by fitting the SWCC with the Gardner model as shown in Figure 4. It is found that the [[alpha].sub.g] = 1.6 x [10.sup.-]2 and the MAE was 1.78 x [10.sup.-]2.

Using the Gardner model, the nonlinear relationship between the matric suction and unsaturated hydraulic conductivity can be obtained as shown in Figure 5. In this example, the total simulation time is one hour. Figure 6 depicted the comparison of the solution with the numerical solution using the Skeel and Berzins' method (Skeel, Berzins 1990). It is found that the computed results agreed well with those obtained from the Skeel and Berzins' method and the MAE was 2.74 x [10.sup.-4] m.

2.3. Verification example 3

The governing equation for unsaturated flow in response to infiltration at ground surface is known as the Richards equation. For one-dimension, it can be written as Eq. (9). Iverson assumed that the Richards equation as a linear diffusion equation and derived the analytical solution (Iverson 2000). Baum et al. revised the analytical solution and implemented in TRIGRS (Baum et al. 2002). The analytical solution proposed by Baum et al. can be written as follows:

[mathematical expression not reproducible], (49)

where Z = z/cos[alpha], d is the depth of steady-state water table measured in z direction, [beta] = [lambda]cos[alpha] and [lambda] = cos[delta] - [[[I.sub.z] / [K.sub.z]].sub.LT], where [I.sub.z] is the long term rainfall flux at ground surface. [K.sub.z] is the hydraulic conductivity, [I.sub.nz] is the surface flux for the nth time interval, n is the total number of time intervals, [D.sub.1] is related to the hydraulic diffusivity [D.sub.0], H(t - [t.sub.n+l]) is the Heaviside step function. The ierfc is the complementary error function and can be expressed as follows:

ierfc([eta]) = 1/[square root of [pi]]exp(-[[eta].sup.2]) - [eta]erfc([eta]). (50) VTC

We conducted a comparison example with the consideration of unsaturated flow in response to infiltration at ground surface. The input parameters were listed in Table 2. The numerical results of this example were compared to those obtained from the analytical solution of equation. In a scenario set in this study, a soil layer thickness is assumed to be 100 cm and a water table at the junction between the soil base and rock bed. Figure 7 and Figure 8 depict that the results from the proposed model and the analytical solution. It is interested that the analytical solution from the simplified linear diffusion equation cannot obtain the appropriate result when the value of the hydraulic diffusivity is large as revealed in Figure 8. The results from this example demonstrate that the analytical solution of equation can only be applied on the problems within certain condition and should be used with caution.

3. Application examples

3.1. Application example 1

The first application example under investigation is a one-dimensional, transient, unsaturated groundwater flow problem in a homogeneous soil. The thickness of the soil L is 150 cm and the slope angle is 31 degrees. Since the Richards equation is an initial value problem, an initial condition must be provided to solve the equation. The initial pressure head is given as follows:

h (z, t = 0) = (-z - [z.sub.w]) cos [alpha], (51)

where [z.sub.w] represents the height of water table. The pore pressure head distribution in unsaturated zones is congruent with the water table function. In this example, we considered [z.sub.w] = 50 cm. The total simulation time is five hours. Considering that a pool of water at the ground surface is maintained holding the infiltration rate, the boundary condition can be defined according to Darcy's law q = -[K.sub.s][nabla]H, and assumes that infiltration rate facilitates complete infiltration and q = [R.sub.f] as follows:

[R.sub.f] = -[K.sub.s] [partial derivative]H/[partial derivative]z. (52)

The definition of [K.sub.s] and H are the same as defined in the previous section. Using Eqs. (4) and (5), the above equation can be written as:

[R.sub.f] = -[K.sub.s] ([partial derivative]h/[partial derivative]z + cos [alpha]). (53)

The above equation can then be rewritten as the following equation which is also known as the Neumann boundary condition:

[partial derivative]h/[partial derivative]z = - cos [alpha] - [R.sub.f]/[K.sub.s]. (54)

If the infiltration rate is greater than saturated hydraulic conductivity, the maximal soil infiltration equals the saturated hydraulic conductivity. The bottom boundary may be assumed to be a no-flow boundary, and is inferred as follows:

0 = -[K.sub.s] [partial derivative]H/[partial derivative]z. (55)

The above equation can be written as:

0 = -[K.sub.s] ([partial derivative]h/[partial derivative]z + cos [alpha]). (56)

Then, we have

[partial derivative]h/[partial derivative]z = -cos [alpha]. (57)

To solve the Richards equation, the boundary conditions must be considered. We adopted [R.sub.f] = 5 x [10.sup.-5] cm/hr. The soil was assumed to be sand, which has saturated hydraulic conductivity ([K.sub.s]) of 9 x [10.sup.-2] cm/hr. The saturated water content ([[theta].sub.s]) and residual water content ([[theta].sub.r]) are 0.43 and 0.08, respectively. The [[alpha].sub.g] of the Gardner model is 1 x [10.sup.-2]. In this example, the total simulation time is five hours. The effective cohesiveness, effective friction angle and unit weight of soil are 4.6 kPa, 30 degrees and 21.5 kN/[m.sup.3], respectively. Figure 9 and Figure 10 indicate the computed results of pressure head and FS for application example 1.

3.2. Application example 2

The second application example under investigation is a one-dimensional, transient, unsaturated groundwater flow problem in two-layered soil. The thickness of the soil L is 150 cm. In the case of two layered soils, layer A and layer B have thickness of 130 cm and 20 cm, respectively, as shown in Figure 11. Table 3 shows the soil parameters for two soil layers. In this example, layer B is more permeable than layer A. The hillslope has a slope angle of 31 degrees. The initial pressure head can be described as Eq. (53). The boundary conditions of the top and bottom boundaries are showed as Eq. (56) and Eq. (59). In this example, we considered [z.sub.w] = 80 cm and [R.sub.f] = 5 x [10.sup.-5] cm/hr. The total simulation time is five hours. The computed pressure heads versus elapsed time for this example are depicted in Figure 12. Figure 12 demonstrates that the fluctuation of pore water pressure in unsaturated layered soil dominates slope stability behavior. Besides, the pressure head at the interface increases quickly and the lowest FS may occur at the interface between two consecutive soil layers during a rainfall even if the top soil layer has the lower hydraulic conductivity. Figure 13 shows the computed FS with respect to time. It shows that the lowest FS occurs at the interface between two consecutive soil layers during a rainfall event. Figure 14 indicates the results of the FS at the depth of 135 cm for this example. It is found that the FS is strongly affected by the rainfall with respect to time. Accordingly, we may conclude that the variation of the pressure head caused by the rainfall event is strongly associated with soil layers and the interface between two consecutive soil layers may play a crucial role for the slope stability.

4. Discussions

In this paper, we carried out three verification examples and two application examples. From the results of verification examples, it is found that the computed numerical results agreed well with the analytical transient solution and those from the Skeel and Berzins' method. Comparison between the proposed numerical model and TRIGRS model (a simplified analytical-based model by assuming the Richards equation as a linear diffusion equation; Iverson 2000; Baum et al. 2002) was also conducted. We conducted a comprehensive study in this example with the consideration of several scenarios. It is found that the proposed model in this study can obtain more appropriate results when the value of the hydraulic diffusivity is larger than those from the TRIGRS model.

In the application examples, we further adopted the proposed model to study the slope stability of a landslide with the consideration of the infiltration from the ground surface for homogenous and layered soils. Since the numerical modeling of rainfall-induced shallow landslides in unsaturated layered soil using the variably saturated flow equation is still hardly found, this study demonstrates the advantages of using the proposed model to study the slope stability of landslides for layered soils. From the results of application examples, we also found that the fluctuation of pore water pressure in unsaturated layered soil dominates slope stability of landslides. Accordingly, it can be concluded that the proposed model in this study may have great potential to deal with problems of rainfall-induced shallow landslides in unsaturated layered soil.

The nature of the Richards equation is highly nonlinear and cannot directly provide an analytical solution. The SWCC is a major factor that influences the nonlinear physical relationship. To derive the proposed numerical model of rainfall-induced shallow landslides in unsaturated layered soil, we adopted the Gardner model to formulate the linearized Richards equation. It should be noted that the Gardner model is based on the one-parameter equation with an indication of the rate of desaturation of a soil. Parameters used in mathematical models for the SWCC include fixed points and two or more fitting constants that are usually selected to capture the general shape of the SWCC. Nevertheless, the Gardner model has only one parameter which may have limitations for fitting the SWCC with complicated shape.

Conclusions

In this study, a numerical model of unsaturated layered soil for rainfall-induced shallow landslides was developed. The infinite slope stability analysis coupled with the hydrological model with the consideration of the SWCC proposed by Gardner was developed to model the shallow landslides triggered by rainfall. The fundamental concepts and the construct of the proposed method are clearly addressed. Findings from this study are drawn as follows:

1. The FDM for the numerical discretization of the variably saturated flow equation based on the Gardner model to deal with groundwater flow in unsaturated layered soil is derived. The proposed method can be used to solve slope stability problems with the appearance of layered soil.

2. The validity of the model is established for a number of test problems by comparing numerical results with the analytical solutions. It is found that the proposed method can be used to model groundwater flow in unsaturated layered soil with high accuracy.

3. The results obtained from this study demonstrate that the slope stability of landslides is strongly dependent to the hydraulic conductivity. It is found that the variation of pore water pressure in unsaturated layered influences the stability of a slope. Besides, the pressure head at the interface increases quickly and the lowest FS may occur at the interface between two consecutive soil layers during a rainfall even if the top soil layer is less permeable.

https://doi.org/10.3846/16486897.2017.1326925

Acknowledgements

This study was partially supported by the Central Geological Survey, Ministry of Economic Affairs, R.O.C. The authors would like to thank the Central Geological Survey for the generous financial support.

Disclosure statement

The authors declare that they don't have any competing financial, professional, or personal interests from any other parties.

References

Alvioli, M.; Baum, R. L. 2016. Parallelization of the TRIGRS model for rainfall-induced landslides using the message passing interface, Environmental Modelling & Software 81: 122-135.

Aylor, D. E.; Parlange, J. Y. 1973. Vertical infiltration into a layered soil, Soil Science Society of America Journal 37(5): 673-676. https://doi.org/10.2136/sssaj1973.03615995003700050015x

Baum, R. L.; Godt, J. W.; Savage, W. Z. 2010. Estimating the timing and location of shallow rainfall-induced landslides using a model for transient, unsaturated infiltration, Journal of Geophysical Research 115(F03013). https://doi.org/10.1029/2009JF001321

Baum, R. L.; Savage, W. Z.; Godt, J. W. 2002. TRIGRS--A FORTRAN program for transient rainfall infiltration and grid-based regional slope-stability analysis, U.S. Geological Survey Open-File Report 02-0424.

Bordoni, M.; Meisina, C.; Valentino, R.; Bittelli, M.; Chersich, S. 2015. Site-specific to local-scale shallow landslides triggering zones assessment using TRIGRS, Natural Hazards and Earth System Science 15(5): 1025-1050. https://doi.org/10.5194/nhess-15-1025-2015

Corradini, C.; Melone, F.; Smith, R. E. 2000. Modeling local infiltration for a two-layered soil under complex rainfall patterns, Journal of Hydrology 237: 58-73. https://doi.org/10.1016/S0022-1694(00)00298-5

Crosta, G. B.; Frattini, P. 2003. Distributed modeling of shallow landslides triggered by intense rainfall, Natural Hazards and Earth System Science 3: 81-93. https://doi.org/10.5194/nhess-3-81-2003

Csete, M.; Buzasi, A. 2016. Climate-oriented assessment of main street design and development in Budapest, Journal of Environmental Engineering and Landscape Management 24(4): 258-268. https://doi.org/10.3846/16486897.2016.1185431

D'Odorico, P.; Fagherazzi, S.; Rigon, R. 2005. Potential for landsliding: dependence on hyetograph characteristics, Journal of Geophysical Research 110(F01007). https://doi.org/10.1029/2004JF000127

Das, B. M. 2010. Principles of geotechnical engineering. 7th ed. Cengage Learning. USA.

Duffy, D. G. 1994. Transform methods for solving partial differential equations. CRC Press, Boca Raton, Fla.

Fok, Y. S. 1970. One-dimensional infiltration into layered soils, Journal of the Irrigation and Drainage Division 96(2): 121-129.

Frattini, P.; Crosta, G. B.; Fusi, N.; Negro, P. D. 2004. Shallow landslides in pyroclastic soil: a distributed modeling approach for hazard assessment, Engineering Geology 73: 277295. https://doi.org/10.1016Zj.enggeo.2004.01.009

Fredlund, D. G.; Morgenstern, N. R. 1977. Stress state variables for unsaturated soils, Journal of the Geotechnical Engineering Division 103: 447-466.

Fredlund, D. G.; Morgenstern, N. R.; Widger, R. A. 1978. The shear strength of unsaturated soils, Canadian Geotechnical Journal 15(3): 313-321. https://doi.org/10.1139/t78-029

Fredlund, D. G.; Xing, A.; Huang, S. 1994. Predicting the permeability function for unsaturated soils using the soil-water characteristic curve, Canadian Geotechnical Journal 31(4): 533-546. https://doi.org/10.1139/t94-062

Gardner, W. R. 1958. Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table, Soil science 85(4): 228-232. https://doi.org/10.1097/00010694-195804000-00006

Gariano, S. L.; Guzzetti, F. 2016. Landslides in a changing climate, Earth-Science Reviews 162: 227-252. https://doi.org/10.1016/j.earscirev.2016.08.011

Green, W. H.; Ampt, G. A. 1911. Studies on soil physics I. The flow of air and water through soils, Journal of Agricultural Science 4: 1-24.

Hachum, A. Y.; Alfaro, J. F. 1980. Rain infiltration into layered soils: prediction, Journal of the Irrigation and Drainage Division. American Society of Civil Engineers 106: 311-319.

Hanks, R.; Bowers, S. 1962. Numerical solution of the moisture flow equation for infiltration into layered soils, Soil Science Society of America Journal 26(6): 530-534. https://doi.org/10.2136/sssaj1962.03615995002600060007x

Iverson, R. M. 2000. Landslide triggering by rain infiltration, Water Resources Research 36(7): 1897-1910. https://doi.org/10.1029/2000WR900090

Keim, R. F.; Skaugset, A. E. 2003. Modelling effects of forest canopies on slope stability, Hydrological Processes 17: 1457-1467. https://doi.org/10.1002/hyp.5121

Ku, C. Y.; Tsai, Y. H. 2013. Solving nonlinear problems with singular initial conditions using a perturbed scalar homotopy method. International Journal of Nonlinear Sciences and Numerical Simulation 14(6): 367-375. https://doi.org/10.1515/ijnsns-2013-0029

Lan, H. X.; Lee, C. F.; Zhou, C. H.; Martin, C. D. 2005. Dynamic characteristics analysis of shallow landslides in response to rainfall event using GIS, Environmental Geology 47: 254-267. https://doi.org/10.1007/s00254-004-1151-8

Liu, C. N.; Wu, C. C. 2008. Mapping susceptibility of rainfall-triggered shallow landslides using a probabilistic approach, Environmental Geology 55(4): 907-915. https://doi.org/10.1007/s00254-007-1042-x

Liu, C. Y.; Ku, C. Y.; Huang, C. C.; Lin, D. G.; Yeih, W. 2015. Numerical solutions of groundwater flow in unsaturated layered soil with extreme physical property contrasts, International Journal of Nonlinear Sciences and Numerical Simulation 16: 325-335. https://doi.org/10.1515/ijnsns-2015-0060

Lu, N.; Godt, J. 2008. Infinite slope stability under steady unsaturated seepage conditions, Water Resources Research 44(W11404). https://doi.org/10.1029/2008WR006976

Lu, N.; Likos, W. J. 2004. Unsaturated soil mechanics. John Wiley and Sons, Inc.

Moldrup, P.; Rolston, D. E.; Hansen, L. A. 1989. Rapid and numerically stable simulation of one dimensional, transient water flow in unsaturated, layered soils, Soil Science 1483(3): 219-226. https://doi.org/10.1097/00010694-198909000-00009

Montgomery, D. R.; Dietrich, W. E. 1994. A physically based model for the topographic control on shallow landslide, Water Resources Research 30(4): 1153-1171. https://doi.org/10.1029/93WR02979

Pereira, P.; Oliva, M.; Baltrenaite, E. 2010. Modelling extreme precipitation in hazardous mountainous areas, Journal of Environmental Engineering and Landscape Management 18(4): 329-342. https://doi.org/10.3846/jeelm.2010.38

Romano, N.; Brunoneb, B.; Santini, A. 1988. Numerical analysis of one-dimensional unsaturated flow in layered soils, Advances in Water Resources 21(4): 315-324. https://doi.org/10.1016/S0309-1708(96)00059-0

Samani, Z.; Cheraghi, A.; Willardson, L. 1989. Water movement in horizontally layered soils, Journal of irrigation and drainage engineering 115(3): 449-456. https://doi.org/10.1061/(ASCE)0733-9437(1989)115:3(449)

Sarkar, K.; Singh, T. N.; Verma, A. K. 2012. A numerical simulation of landslide-prone slope in Himalayan region--a case study, Arabian Journal of Geosciences 5(1): 73-81. https://doi.org/10.1007/s12517-010-0148-8

Schiliro, L.; Esposito, C.; Mugnozza, G. S. 2015. Evaluation of shallow landslide-triggering scenarios through a physically based approach: an example of application in the southern Messina area (northeastern Sicily, Italy), Natural Hazards and Earth System Science 15(9): 2091-2109. https://doi.org/10.5194/nhess-15-2091-2015

Singh, T. N.; Verma, A. K.; Sarkar, K. 2010. Static and dynamic analysis of a landslide, Geomatics, Natural Hazards and Risk 1(4): 323-328. https://doi.org/10.1080/19475705.2010.521354

Skeel, R. D.; Berzins, M. 1990. A method for the spatial discretization of parabolic equations in one space variable, SIAM Journal on Scientific and Statistical Computing 11(1): 1-32. https://doi.org/10.1137/0911001

Srivastava, R.; Yeh, T. C. 1991. Analytical solutions for one-dimensional, transient infiltration toward the water table in homogeneous and layered soils, Water Resources Research 27(5): 753-762. https://doi.org/10.1029/90WR02772

Tothova, I.; Igaz, D.; Antal, J. 2007. Soil water movement modelling in hapllicc luvisols (HMm) and albi hapllic luvisols (HMl) under Slovak climatic conditions, Journal of Environmental Engineering and Landscape Management 15(2): 69-75.

Tracy, F. T. 2006. Clean two- and three-dimensional analytical solutions of Richards' equation for testing numerical solvers, Water Resources Research 42: W08503. https://doi.org/10.1029/2005WR004638

Tracy, F. T. 2011. Analytical and numerical solutions of Richards' equation with discussions on relative hydraulic conductivity. INTECH Open Access Publisher, Janeza Trdine 9, 51000 Rijeka, Croatia.

Tsai, T. L.; Chen, H. F. 2010. Effects of degree of saturation on shallow landslides triggered by rainfall, Environmental Earth Sciences 59(6): 1285-1295. https://doi.org/10.1007/s12665-009-0116-3

Tsai, T. L.; Yang, J. C. 2006. Modeling of rainfall-triggered shallow landslide, Environmental Geology 50(4): 525-534. https://doi.org/10.1007/s00254-006-0229-x

Van Genuchten, M. T. 1980. A closed form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Science Society of American Journal 44(5): 892-898. https://doi.org/10.2136/sssaj1980.03615995004400050002x

Vanapalli, S. K.; Fredlund, D. G.; Pufahl, D. E.; Clifton, A. W. 1996. Model for the prediction of shear strength with respect to soil suction, Canadian Geotechnical Journal 33(3): 379-392. https://doi.org/10.1139/t96-060

Wang, C.; Esaki, T.; Xie, M.; Qiu, C. 2006. Landslide and debris-flow hazard analysis and prediction using GIS in Minamata-Hougawachi area, Japan, Environmental Geology 51: 91-102. https://doi.org/10.1007/s00254-006-0307-0

Whisler, F.; Klute, A. 1966. Analysis of infiltration into stratified soil columns, in Proc Wageningen Symp, Proc. IAHS Symposium on Water in the Unsaturated Flow. Wageningen, The Netherlands, 451-470.

Yeh, H. F.; Lee, C. H. 2013. Soil water balance model for precipitation-induced shallow landslides, Environmental Earth Sciences 70(6): 2691-2701. https://doi.org/10.1007/s12665-013-2326-y

Yoshimatsu, H.; Abe, S. 2006. A review of landslide hazards in Japan and assessment of their susceptibility using an analytical hierarchic process (AHP) method, Landslides 3: 149-158. https://doi.org/10.1007/s10346-005-0031-y

Chih-Yu LIU. Doctoral Candidate, Department of Harbor and River Engineering, National Taiwan Ocean University, Taiwan, 2 peer-reviewed articles in scientific journals, 14 conference papers.

Cheng-Yu KU. PhD, Professor, Department of Harbor and River Engineering, National Taiwan Ocean University, Taiwan, 58 peer-reviewed articles in scientific journals, 4 book chapters, 4254 citations, H-index 18 (Google Scholar).

Jing-En XIAO. Doctoral Candidate, Department of Harbor and River Engineering, National Taiwan Ocean University, Taiwan, 1 peer-reviewed articles in scientific journals, 13 conference papers.

Chi-Chao HUANG. Doctoral Candidate, Department of Resources Engineering, National Cheng Kung University, Taiwan, Chief Secretary, Central Geological Survey, Ministry of Economic Affairs, Taiwan.

Shih-Meng HSU. PhD, Assistant Professor, Department of Harbor and River Engineering, National Taiwan Ocean University, Taiwan, 32 peer-reviewed articles in scientific journals, 3 book chapters, 254 citations, H-index 8 (Google Scholar).

Chih-Yu LIU (a), Cheng-Yu KU (a), Jing-En XIAO (a), Chi-Chao HUANG (b), Shih-Meng HSU (a)

(a) Department of Harbor and River Engineering, National Taiwan Ocean University, Keelung, Taiwan

(b) Central Geological Survey, Ministry of Economic Affairs, Taiwan

Submitted 05 Jan. 2017; accepted 022 May 2017

Corresponding author: Cheng-Yu Ku

E-mail: chkst26@maiSntou.edu.tw

Caption: Fig. 1. The definition of the coordinate system in this study

Caption: Fig. 2. The configuration of water flow through layered soils in the FDM formulation

Caption: Fig. 3. Result comparison of verification example 1 with the exact solution

Caption: Fig. 4. The fitting curve of the SWCC using the Gardner model

Caption: Fig. 5. The relationship between matric suction and unsaturated hydraulic conductivity for the Gardner model

Caption: Fig. 6. Result comparison of verification example 2 with the numerical solution

Caption: Fig. 7. Comparison of the results for [D.sub.0] = 10 x [K.sub.s] (left: this study, right: TRIGRS model)

Caption: Fig. 8. Comparison of the results for [D.sub.0] = 100 x [K.sub.s] (left: this study, right: TRIGRS model)

Caption: Fig. 9. Computed pressure head of the soil profile for application example 1

Caption: Fig. 10. Computed FS of the soil profile for application example 1

Caption: Fig. 11. Schematic illustration of application example 2

Caption: Fig. 12. Computed pressure head of the soil profile for application example 2

Caption: Fig. 13. Computed FS of the soil profile for application example 2

Caption: Fig. 14. The computed results of FS at z = 135 cm for application example 2

Table 1. Parameters used in the verification example 1 saturated hydraulic conductivity (m/s) 9 x [10.sup.-5] [[alpha].sub.g] of the Gardner model [10.sup.-4] saturated water content 0.50 residual water content 0.11 Table 2. Parameters used in the verification example 3 This study Baum et al. (2002) Slope angle (degree) 30 Thickness (m) 1 Inflitration rate (m/hr) 0.2 Inflitration time (hr) 5 Saturated hydraulic 2.5 conductivity (m/hr) [[theta].sub.s] = 0.43 [D.sub.0] = (10~100) x Other parameters [[theta].sub.r] = 0.08 [K.sub.s] [[alpha].sub.z] = 0.2 (m/hr) Table 3. Parameters used in the application example 2 Soil Layer A Layer B Hydraulic conductivity 1 x [10.sup.-2] 1 x [10.sup.-1] [K.sub.s] (cm/hr) Saturated water content 0.43 0.50 [[theta].sub.2] Residual water content 0.08 0.11 [[theta].sub.r] [[alpha].sub.g] of 1 x [10.sup.-2] 1 x [10.sup.-2] the Gardner model

Please Note: Illustration(s) are not available due to copyright restrictions.