Application of the Chebyshev spectral method to the simulation of groundwater flow and rainfall-induced landslides

Journal Pre-proof

Application of the Chebyshev spectral method to the simulation of groundwater flow and rainfall-induced landslides L.Z. Wu , S.R. Zhu , P.J. Peng PII: DOI: Reference:

S0307-904X(19)30726-7 https://doi.org/10.1016/j.apm.2019.11.043 APM 13173

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

26 July 2019 20 November 2019 26 November 2019

Please cite this article as: L.Z. Wu , S.R. Zhu , P.J. Peng , Application of the Chebyshev spectral method to the simulation of groundwater flow and rainfall-induced landslides, Applied Mathematical Modelling (2019), doi: https://doi.org/10.1016/j.apm.2019.11.043

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Inc.

Highlights • CSM was developed to model groundwater flow in unsaturated soils. • CSM has higher precision than FDM and FEM with fewer grid points. • CSM successfully solves rainfall infiltration problems for an unsaturated slope.

1

Application of the Chebyshev spectral method to the simulation of groundwater flow and rainfall-induced landslides L. Z. Wu a*, S. R. Zhu a, P. J. Peng a, b a) State Key Laboratory of Geohazard Prevention and Geoenvironment Protection, Chengdu University of Technology, 610059, Chengdu, Sichuan, P. R. China b) College of Geological Engineering and Geomatics, Chang'an University, Xi'an, 710054, Shaanxi, China

* Corresponding author E-mail addresses: [email protected] (L. Z. Wu), [email protected] (S. R. Zhu), [email protected] (J. B. Peng)

Resubmitted on November 20, 2019

2

Abstract In this paper, the potential of the Chebyshev spectral method (CSM) is examined for the numerical solution of time-dependent variably saturated Darcian flow problems described by the Richards equation (RE). Generally, the computational efficiency of the traditional finite difference method (FDM) is not high because it requires a high mesh density to improve accuracy. Therefore, the Chebyshev spectral method (CSM) was extended to model groundwater transient flow in unsaturated soils. Furthermore, the Gardner model (soil-water characteristic curve) was adopted to linearize the RE for an inclined slope. The test results demonstrated that the accuracy, efficiency, and robustness of the CSM are better than those of the traditional FDM. A comparison of the analytical solutions of each method showed that the computational accuracy ( L2  h*  ) of the CSM was less affected by the grid size than that of the FDM. Simultaneously, the CSM was not sensitive to the initial conditions. The results demonstrated that the numerical accuracy of the CSM was stable in the order of 10-6 to 10-7, whereas that of the FDM was in the order of 10-3 to 10-6; it achieved higher accuracy with fewer grid points. The CSM was successfully applied to solve the one-dimensional transient flow problem on unsaturated soil slopes. The numerical results indicated that the proposed method solved the transient flow problem related to rainfall-induced landslides with high accuracy. This proposed approach can be developed to analyze rainfall-induced landslides. Keywords: Chebyshev spectral method; Richards equation; slope stability; accuracy; grid size. 3

1. Introduction Global warming is increasing the number of heavy rainfall events, which is leading to many geological hazards [1-5]. Rainfall infiltration affects slope stability, thus it is of practical significance to develop a rainfall infiltration equation and evaluate the unsaturated soil slope performance [6-10]. The Richards equation (RE) [11] is commonly used to describe rainfall infiltration in unsaturated soils; however, it is difficult to obtain analytical solutions [12-14] because the high nonlinearity of the RE is related to the complex hydraulic properties of unsaturated soils [15-17]. Hence, the numerical solution of the RE has been extensively studied to examine groundwater flow [18-23]. A number of advanced methods have emerged to solve variably saturated flow problems, such as isogeometric analysis [24], the finite volume method [25], the meshless method [26, 27], and some high-order methods [28, 29]. For example, a procedure

for

variably

saturated

flow

in

porous

media

with

complex

three-dimensional (3D) geometries was developed using a 3D finite volume unstructured mesh modular framework [30]. A coordinate transform of nonorthogonal grids was used to solve saturated-unsaturated flow problems in complex terrain [31]. An approximation solution was obtained for a two-dimensional (2D) time-dependent unsaturated flow equation using second-order accurate cell-centered finite-volume discretization on unstructured grids [32]. Recently, transient flow in unsaturated soils was solved using a novel collocation meshless method [33]. A novel generalized finite difference scheme was developed to solve the RE based on nonrectangular structured 4

grids [34]. An adaptive higher-order space-time discontinuous Galerkin method was presented for solving the RE to optimize accuracy and efficiency [29]. Although the use of improved methods for solving unsaturated flow is successful from different aspects, there is a growing need to develop an advanced computing method to improve numerical accuracy and robustness [35, 36]. Consequently, many advanced numerical methods have been proposed, such as the Chebyshev spectral method (CSM) [37], Trefftz method [26], meshless method based on fundamental solutions [27], and diffuse element method [38]. Among these methods, the CSM can solve partial differential equations effectively. The original spectral method was developed to improve the accuracy of the numerical simulation of sound wave propagation [39, 40]. Numerical approaches such as the finite element method [41] and finite difference method (FDM) [42] are widely used for the discrete RE to obtain linear equations. A traditional numerical scheme, such as the FDM, is constructed based on the Crank– Nicolson FDM in time [43]. However, the computational efficiency is not high because it requires a high mesh density to improve accuracy [44]. Therefore, the CSM will be helpful in solving the RE to obtain higher accuracy and better stability [45]. Currently, no published studies have used the CSM to solve transient flow in unsaturated soils and there are no findings on the CSM’s application to the analysis of rainfall-induced landslides. In this study, a new numerical method is developed for solving the RE. The proposed method with the CSM and other analytical solutions are compared, and the 5

robustness and accuracy of the proposed method are validated. The effect of the number of nodes or grid sizes and initial conditions on the results of the CSM and traditional FDM are evaluated. Additionally, solutions for transient flow in an unsaturated slope are also obtained to demonstrate the proposed method.

2. Numerical procedures 2.1. Governing equation and finite difference method The 2D generalized RE for groundwater flow in unsaturated soils is expressed as [11]   H  Kx  h x  x

H       Kz h z  z 

H  ,   C h t 

(1)

where K z (h) and K x (h) are the permeability coefficients along the vertical direction and lateral direction in unsaturated soils, respectively; H is the total head; h is the pressure head; and C h  is the specific moisture capacity, which can be defined by

C  h    / h . The original RE is only horizontal. To study the groundwater flow of the unsaturated slope (Fig. 1), the RE needs to be rotated. Total head H can be written as H  hE ,

(2)

and elevation head E can be described as

E  x sin    z cos   ,

(3)

where  is the slope angle. By substituting Eqs. (2) and (3) into Eq. (1), Eq. (1) can be re-expressed as

6

  h  h     h  K x  h    sin        K z  h    cos       C  h  .  x  t  x   z   z 

(4)

To solve Eq. (4), the permeability coefficient and moisture content in unsaturated soils are assumed to follow the exponential form: K  h  Kse

g h

,

(5)

  h   r   s  r  e

g h

,

(6)

where  s and  r are the saturated and residual water content, respectively; K s is the saturated permeability coefficient of soils; and  g is a pore-size distribution parameter [46]. According to Iverson’s model [47], the modified RE only considers flow in the vertical direction, which can be given by

  h    K z  h    cos       .  z   z   t

(7)

The RE for an inclined slope can be obtained by substituting Eqs. (5) and (6) and a new parameter h* into Eq. (7) as follows [42]:

 2 h* h* h*   cos   c ,   g z 2 z t  h  where c   g  s  r  / K s , h*  e  e g

g hd

(8)

, and hd is the pressure head for dry

soil. The FDM is adopted to obtain an approximate solution to Eq. (8). Simultaneously, the Crank–Nicolson implicit difference format can be expressed as 1  hi*1n  2hi*n  hi*1n hi*1n 1  2hi*n 1  hi*1n 1     z 2  2 2  ag cos( )  hi*1n  hi*1n hi*1n 1  hi*1n 1   hi*n  hi*n 1      c  2z  2 2 t   .

(9)

where i, z, t denote the nodal point, discrete interval along the z-axis, and time interval, respectively, and n denotes the time step. Moreover, the simplified matrix of 7

Eq. (9) can be expressed as Ah*  b , where A is a h* is a

 N  1 1

 N  1 1

vector, and b is a

 N  1   N  1

square matrix,

vector, where N+1 is the number

of nodes. The linear equations can be solved using the triangular factorization algorithm.

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

z j  cos  j / N  , j  0,1,

N.

(10)

The first-order derivative of interpolation function p  z  can be written as

p  z   DN h , where z

(11)

is represented as the N+1 dimensional vector

represents the N+1 dimensional vector head, and D N is an

 N  1   N  1

 h0 , h1 ,

, hN 

T

 z0 , z1 ,

, zN  , h T

composed of the pressure

Chebyshev differentiation matrix.

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

2N 2 1 ,  DN 00  6

 DN  jj 

z j

2 1  z 2j 

 DN NN

2N 2 1  , 6

j  1,

,

8

, N 1 ,

(12a)

(12b)

 DN ij

c  1  i , c j zi  z j i j

(i  j )i, j  0, , N ,

(12c)

where  D N ij represents the element in row i+1 and column j+1 of Chebyshev differentiation matrix D N for the first derivative and

i  0, N . i  1, , N  1

2, ci   1,

(13)

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

  DN , z

(14.a)

2  D2N , 2 z

(14.b)

3  D3N . z 3

(14.c)

Because D N is obtained in the interval [-1, 1], it is necessary to scale D N to other intervals, that is,  2  DN z L

(15.a)

n 2    DN  . n z L  n

(15.b)

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

1 2 * g h* DN h  cos    DN h*  . c c t

(16)

As shown in Eq. (16), Eq. (8) has been transformed into an ordinary differential equation (ODE) by introducing Chebyshev differentiation matrices. In this study, the 9

variable-step four- and fifth-order Runge–Kutta method is used to solve Eq. (16), which is the ODE solver in MATLAB. Meanwhile, the 2D linear RE can also be solved using the same method. To evaluate the performance of the proposed method, an L2 norm is used to indicate the error of solution h* :

L2  h*   h*,ana ( z, t )  h* ( z, t ) , 2

(17)

where h*,ana ( z, t ) and h* ( z, t ) are one-dimensional (1D) analytical and numerical solutions, respectively. Simultaneously, the numerical and analytical solutions with spatially and temporally grids can produce an absolute error defined as

errA  h*   h*,ana ( z, t )  h* ( z, t ) .

(18)

3. Verification of the methods 3.1. Test 1 Test 1 describes transient groundwater flow in homogeneous unsaturated soils. The mathematical model is shown in Fig. 2. The soil thickness (L) is assumed to be 10 m. The model parameters are  g  1104 ,  s  0.50 , r  0.11 , and K s =9×10-5 m/h [22]. The governing equation is described as

 h  K z  h  h  C h .  Kz  h   z  z  z t

(19)

When water infiltrates into a soil mass, ponding on the ground maintains the pressure head at zero [48]. As shown in Fig. 2, the upper and lower boundary conditions can be described as [22]

h  z  0, t   hd , 10

(20)

h  z  L, t   0 .

(21)

The normalized analytical solution for the pressure head can be written as [49]

h*  z, t   ht *  z, t   hs*  z  ,

(22)

where hs*  z  is the analytical solution for steady-state flow,



hs*  z   1  e

 g hd

  1  e  / 1  e  ,  g z

 g L

(23)

and ht*  z, t  is the analytical solution for transient flow,

ht  z, t   *



 g hd

2 1 e

Lc



e

 g  L  z  /2



  1 m 1

m

 m   m t   sin  m z  e ,   m

(24)



where m  k / L , m   g2 / 4  m2 / c , and t is time. The final analytical solution can be written as h  z, t  

1

g

 h  z, t   e  . *

 g hd

(25)

The effect of z and N on the results of Test 1 was analyzed. The time step was set to 0.01 h. Different numbers of nodes or grid sizes, that is, N /△z = 20 / 0.5 (m), N /△z = 40 / 0.25 (m), and N /△z = 80 / 0.125 (m), were selected to investigate the influence of the different methods on L2  h*  . The total simulation time was 5 h. Simultaneously, three different initial conditions h0 were chosen to assess their influence on the numerical accuracy and stability of the FDM and CSM. As presented in Fig. 3, L2  h*  of the CSM were consistently smaller than that of the traditional FDM. The numerical accuracy increased as the mesh was refined, and the accuracy of the traditional FDM was significantly affected by the initial conditions. However, the proposed method (CSM) was less affected, and the numerical accuracy was stable in the order of 10-6 to 10-7 in this example. The result shows that the CSM 11

was more robust than the FDM. Fig. 4 shows the comparison of the computed pressure head profiles with different grid sizes at t =2 h. Compared with the analytical solution, it can be easily seen that the numerical accuracy of the CSM was less influenced by the grid sizes than that of the FDM (Fig. 4); that is, the CSM achieved better numerical results with fewer mesh nodes. The pressure head profiles were obtained from the computed results over time under h0 = -100,000 m (Fig. 5). The number of nodes or grid sizes was 40 / 0.25 (m). Fig. 5 illustrates that the numerical solutions of the CSM agree well with the analytical solutions over time, whereas the numerical solutions of the FDM have larger errors than the analytical solutions. In Fig. 6, the minimum absolute errors of the proposed CSM and the FDM are approximately 10-11 and 10-6, respectively. Consequently, the extended CSM is characterized by higher accuracy results than traditional the FDM with fewer nodes.

3.2. Test 2 To further verify the effectiveness of the proposed method, 2D transient infiltration was investigated. The geometry and boundary conditions for this example are illustrated in Fig. 7, where L = 1 m and W = 1 m. The soil was assumed to be silt. The parameters of soil were described by Liu [42], and included  g  8 103 , K s  1104 m / s ,  s  0.35 , and r  0.14 . The governing equation is described as

follows:

  h    h  K z  h  h  C h .  Kx  h    Kz h   x  x  z  z  z t

(26)

The initial condition of the model can be written as h  x, z, t  0   hd , where hd 12

was taken as -1 m. The upper boundary condition was set to a special pressure head [9] as described in Eq. (30), and its distribution is shown in Fig. 8. The boundaries of the model are described as follows:

h  0, z, t   hd

(27)

h W , z, t   hd

(28)

h  x,0, t   hd

(29)

   x   g hd  h h  x, L, t   ln  1  e g d sin  e W  





  / g . 

(30)

The normalized analytical solution ht* x, z, t  and hs* x, z, t  for this model can be expressed as [49]

ht  x, z, t   *



 g hd

2 1 e

Lc

 sin   x  e   W 



hs  x, z, t   1  e *

 g hd

 g  L  z  /2



  1 m 1

e

 g  L z 

2

m

 m   mt   sin  m z  e ,  m 

  x  sinh  1 z  sin   ,  W  sinh  1 L 

(31)

(32)

where 1   g2 / 4   / W  . Similarly, the final analytical solution can be 2

obtained using Eqs. (22) and (25). The time step and total simulation time were taken as 0.01 h and 5 h in this example, respectively. The effect of N /△z on the results of the different methods are discussed. In Fig. 9, L2  h*  decreases as the grid sizes decrease, and L2  h* 

of

the CSM are much lower than that of the FDM, which means that the CSM obtained higher numerical precision for the same mesh density in this example. Fig. 10 compares the computed results of the 2D transient infiltration problem using the CSM and analytical solution. Highly accurate numerical solutions in the order of 10−7 were obtained in this example. The results demonstrated that pressure 13

head contours obtained by the CSM agreed well with the analytical solutions, whether the simulation time was 2 h (Fig. 10a) or 5 h (Fig. 10b).

4. Application example This application example simulated transient infiltration in an unsaturated slope using the CSM. The sloped thickness (L) and angle (  ) were assumed to be 10 m and 33°, respectively (Fig. 11). The initial condition was unchanged at h( z, t  0)  10 m. The governing equation is described in Eq. (7). The upper and lower boundary conditions of the soil can be expressed as Eqs. (20) and (21). The exact solution to this example is written as [49]

ht  z, t   *



 g hd

2 1 e

Lc

e

 g cos( ) L  z  /2



  1 m 1

m

 m   mt   sin  m z  e .  m 

(33)

The steady-state solutions hs*  z  of Eq. (8) can be obtained as



hs*  z   1  e

 g hd

1 e

1 e

 g cos   z

 g cos   L

.

(34)

The final analytical solution is obtained as Eq. (25). For this example, we selected sandy soil and silty loam. The experimental data [50, 51] and fitting curve of the two soils are shown in Fig. 12. The fitted parameters, including K s ,  s ,  r , and  g , are listed in Table 1. The total simulation time was 2 h, and N = 40. The process for solving Eq. (8) is similar to that used for Test 1. The pressure head calculated using the CSM was introduced into the infinite slope stability analysis and the factor of safety (FS) could be obtained [33]: FS 

c   z t cos 2     h w cos 2     tan   z t cos    sin   14

,

(35)

where c denotes the effective cohesion,   denotes the effective friction angle,  t denotes the unit weight of soil,  w denotes the unit weight of water, and z denotes the soil thickness. It was assumed that the unit weight of soil was 19 (kN/m3), and the effective cohesiveness and effective friction angle of sandy soil and silty loam were 0 (kPa), 32 (°) and 10 (kPa), 20 (°), respectively [51]. The revised FS for analyzing slope stability can be obtained as FS 

c

z t cos    sin   



tan   h w   h    r /  s   r  tan    tan   z t cos   sin  

(36)

The profiles of the pressure head of two soils over time were compared with the analytical solutions. Fig. 13 shows that the computed results were in good agreement with the analytical solutions. The numerical findings for sandy soil and silty loam indicate that the accuracy of the absolute error can reach the order of 10-5 and 10-6, respectively (Fig. 14). Fig. 15 shows the computed results of FS for different unsaturated soils. Additionally, the range of FS for different soils is different, and FS for silty loam is significantly smaller than that for sandy soil. Table 2 shows variations in FS at different depths over the infiltration time for different soils. It can be seen that the safety factor decreased with increasing depth and infiltration time. The FS for sandy soil was greater than 1.0 in Table 2, and the FS for silty loam was less than 1.0 at z = 9.263m. Consequently, the stability of the slope of sandy soil was better than that of the slope of silty loam. The comparison of the computed results reveals that rainfall-induced landslides are closely related to soil types.

5. Conclusions 15

The CSM, which is based on the discretization of Chebyshev points, was extended to model groundwater flow in porous media. A Chebyshev differentiation matrix is used to construct an n-th order Chebyshev differentiation matrix conveniently and quickly. Then, it is introduced to transform a partial differential equation into an ODE. It is characterized by simple code, reliable operation, and high accuracy. This is the first time that the CSM has been used successfully to solve 1D and 2D transient flow problems in unsaturated soils. Compared with the traditional FDM, the CSM was more efficient, and it achieved higher accuracy with a lower-resolution grid. The CSM was not sensitive to the initial conditions. Additionally, the results demonstrated that the range of numerical accuracy of the CSM was stable in the order of 10-6 to 10-7, whereas that of the FDM was in the order of 10-3 to 10-6. The numerical results obtained by the proposed method were highly consistent with the analytical transient solution. This illustrated that the proposed method had sufficiently high precision to deal with transient flow issues related to rainfall-induced landslides. The CSM introduced into finite slope stability analysis can be applied to analyze rainfall-induced landslides effectively. This is a novel attempt and challenge, particularly in slope stability analysis in the engineering field.

Acknowledgements This work was supported by the Major Program of National Natural Science Foundation of China (no. 41790441), the National Key R&D Program of China (no. 2018YFC1504702), and the National Natural Science Foundation of China (no. 41672282).

Data availability All data generated or analyzed during this study are included within the article. 16

Conflict of interest The authors declare that they have no conflict of interest.

References [1] Kuo, Y. C., Lee, M. A., & Lu, M. M., 2016. Association of Taiwan’s rainfall patterns with large-scale oceanic and atmospheric phenomena. Atmospheric Research, 180, 200–210. [2] Wu, L.Z., Zhou, Y., Sun, P., Shi, J.S., Liu, G.G., Bai, L.Y., 2017. Laboratory characterization of rainfall-induced loess slope failure. Catena, 150: 1–8. [3] Wu, L.Z., Deng, H., Huang, R.Q., Zhang, L.M., Guo, X.G., & Zhou, Y., 2019. Evolution of lakes created by landslide dams and the role of dam erosion: A case study of the Jiajun landslide on the Dadu River, China. Quaternary International. 503A: 41–50. [4] Jun, B. H. 2016, Numerical Simulation of the Topographical Change in Korea Mountain Area by Intense Rainfall and Consequential Debris Flow. Advances in Meteorology, 9363675. [5] Wu, L. Z., Selvadurai, A. P. S., Zhang, L. M., Huang, R. Q., & Huang, J., 2016. Poro-mechanical coupling influences on potential for rainfall-induced shallow landslides in unsaturated soils. Advances in Water Resources, 98, 114–121. [6] Yeh, H. F., & Lee, C. H., 2013. Soil water balance model for precipitation-induced shallow landslides. Environmental Earth Sciences. 70 (6), 2691–2701. [7] Collins, B. D., & Znidarcic, D., 2004. Stability analyses of rainfall induced landslides. Journal of Geotechnical and Geoenvironmental Engineering, 130(4), 17

362–372. [8] Wu, L. Z., Huang, Jinsong, Fan, W., Li, X., 2020. Hydro-mechanical coupling in unsaturated soils covering a non-deformable structure. Computers and Geotechnics, 117: 103287. [9] Ku, C. Y., Liu, C. Y., Su, Y. , Xiao, J. E., & Huang, C. C., 2017. Transient modeling of regional rainfall-triggered shallow landslides. Environmental Earth Sciences, 76(16), 570. [10] Wu, L. Z., Zhang, L. M., Zhou, Y., & Li, B. E., 2017. Analysis of multi-phase coupled seepage and stability in anisotropic slopes under rainfall condition. Environmental Earth Sciences, 76(14): 469. [11] Richards L A. 1931, Capillary Conduction of Liquids Through Porous Mediums. Physics, 1(5), 318–333. [12] Wu, L. Z., Zhang, L.M., 2009. Analytical solution to 1D coupled water infiltration and deformation in unsaturated soils. International Journal for Numerical and Analytical Methods in Geomechanics. 33(6): 773–790. [13] Wu, L.Z., Zhang, L.M., Huang, R.Q., 2012. Analytical solution to 1D coupled water infiltration and deformation in two-layer unsaturated soils. International Journal for Numerical and Analytical Methods in Geomechanics. 36(6): 798– 816. [14] Romano, N., Brunone, B., Santini, A., 1998. Numerical analysis of one-dimensional unsaturated flow in layered soils. Advances in Water Resources, 21(4), 315–324. 18

[15] Srivastava, R., Yeh, T. C. J., 1991. Analytical solutions for one-dimensional, transient infiltration toward the water table in homogeneous layered soils. Water Resour. Res. 27(5), 753–762. [16] Peranić, Josip, Arbanas, Željko, Cuomo, S., & Maček, Matej., 2018. Soil-water characteristic curve of residual soil from a flysch rock mass. Geofluids, 6297819. [17] Genuchten, V., & Th., M., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils1. Soil Science Society of America Journal, 44(5), 892–898. [18] Fredlund, D.G., Xing, A., Huang, S., 1994. Predicting the permeability function for unsaturated soils using the soil-water characteristic curve. Can. Geotech. J. 31(4), 533–546. [19] Wu LZ, Liu GG, Wang LC, Zhang LM, Li BE, Li B., 2016. Numerical analysis of 1D coupled infiltration and deformation in layered unsaturated porous medium. Environmental Earth Sciences. 75(9), 761. [20] Ku, C. Y., Liu, C. Y., Xiao, J.E., & Yeih, W., 2017. Transient modeling of flow in unsaturated soils using a novel collocation meshless method. Water, 9(12), 954. [21] Zeng, J., Zha, Y., & Yang, J., 2018. Switching the Richards’ equation for modeling soil water movement under unfavorable conditions. Journal of Hydrology, 563: 942–949. [22] Zhu, S.R., Wu, L.Z., Shen, Z.H., & Huang, R.Q., 2019. An improved iteration method for the numerical solution of groundwater flow in unsaturated soils. Computers and Geotechnics. 2019, 114, 103113. 19

[23] Paulus, R., Dewals, B. J., Erpicum, Sébastien, Pirotton, M., & Archambeau, P., 2013. Innovative modelling of 3d unsaturated flow in porous media by coupling independent models for vertical and lateral flows. Journal of Computational and Applied Mathematics, 246, 38-51. [24] Nguyen, M. N., Bui, T. Q., Yu, T., & Hirose, S., 2014. Isogeometric analysis for unsaturated flow problems. Computers and Geotechnics, 62, 257-267. [25] Bevilacqua, I., Canone, D., & Ferraris, S., 2009. Acceleration techniques for the iterative resolution of the Richards equation by the finite volume method. International Journal for Numerical Methods in Biomedical Engineering, 27(8), 1309-1320. [26] Liu, C. S. 2008. A modified collocation Trefftz method for the inverse Cauchy problem of Laplace equation. Engineering Analysis with Boundary Elements, 32(9), 778–785. [27] Xiao, J. E., Ku, C. Y., Liu, C. Y., Fan, C. M., & Yeih, W., 2017. On solving free surface problems in layered soil using the method of fundamental solutions. Engineering Analysis with Boundary Elements, 83, 96–106. [28] Solin, P., & Kuraz, M. 2011. Solving the nonstationary Richards equation with adaptive hp-FEM. Advances in Water Resources, 34(9), 1062–1081. [29] Dolejší, V., Kuraz, M., & Solin, P. 2019. Adaptive Higher-Order Space-Time Discontinuous

Galerkin

Method

for

the

Computer

Simulation

of

Variably-Saturated Porous Media Flows. Applied Mathematical Modelling, 72, 276-305. 20

[30] Mcbride, D., Cross, M., Croft, N., Bennett, C., & Gebhardt, J., 2006. Computational modelling of variably saturated flow in porous media with complex three-dimensional geometries. International Journal for Numerical Methods in Fluids, 50(9), 1085-1117. [31] Deng, B., & Wang, J., 2017. Saturated-unsaturated groundwater modelling using 3D Richards equation with a coordinate transform of nonorthogonal grids. Applied Mathematical Modelling, 50(10), 39-52. [32] Manzini, G., & Ferraris, S. 2004. Mass-conservative finite volume methods on 2-D unstructured grids for the Richards’ equation. Advances in Water Resources, 27(12), 1199–1215. [33] Ku, C. Y., Liu, C. Y., Su, Yan., & Xiao, J. E., 2018. Modeling of transient flow in unsaturated geomaterials for rainfall-induced landslides using a novel spacetime collocation method. Geofluids, 7892789. [34] Chávez-Negrete C, Domínguez-Mota F J, Santana-Quinteros D., 2018. Numerical solution of Richards’ equation of water flow by generalized finite differences. Computers and Geotechnics, 101, 168–175. [35] List, F., & Radu, F. A., 2016. A study on iterative methods for solving Richards’ equation. Computational Geosciences, 20(2), 341–353. [36] Fassino, C., & Manzini, G., 1998. Fast-secant algorithms for the non-linear Richards equation. Communications in Numerical Methods in Engineering, 14(10), 921–930. [37] Heping, Ma, Benyu, & Guo., 1988. The Chebyshev spectral method for 21

burgers-like equations. Journal of Computational Mathematics, 6(1), 48–53. [38] Nayroles, B., Touzot, G., & Villon, P., 1992. Generalizing the finite element method: diffuse approximation and diffuse elements. Computational Mechanics, 10(5), 307-318. [39] Gottlieb, D., Orszag, S. A., & Sod, G. A., 1978. Numerical analysis of spectral methods: theory and application (CBMS-NSF Regional Conference Series in Applied Mathematics). Journal of Applied Mechanics, 45(4), 969. [40] Igel, H. 2010. Wave propagation in three‐dimensional spherical sections by the chebyshev spectral method. Geophysical Journal of the Royal Astronomical Society, 136(3), 559–566. [41] Paniconi, C., & Putti, M., 1994. A comparison of picard and newton iteration in the numerical solution of multidimensional variably saturated flow problems. Water Resources Research, 30(12), 3357–3374. [42] Liu C Y, Ku C Y, Huang C C, et al., 2015. Numerical Solutions for Groundwater Flow in Unsaturated Layered Soil with Extreme Physical Property Contrasts. International Journal of Nonlinear Sciences & Numerical Simulation, 16(7). 325–335. [43] Ravi Kanth, A. S. V., & Garg, N. 2019. An implicit numerical scheme for a class of multi-term time-fractional diffusion equation. The European Physical Journal Plus, 134(6), 312, 1-17. [44] Yang, Z., et al., 2018. Mixed-grid finite-difference methods for wave equation numerical modeling in time-space domain. Lithologic Reservoirs, 30(2), 93-109. 22

[45] Doha, E. H., Bhrawy, A. H., & Ezz-Eldien, S. S., 2011. A Chebyshev spectral method based on operational matrix for initial and boundary value problems of fractional order. Computers & Mathematics with Applications, 62(5), 2364– 2373. [46] 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. [47] Iverson, R. M. 2000. Landslide triggering by rain infiltration. Water Resources Research, 36(7), 1897–1910. [48] Green, W.H., Ampt, G.A. 1911. Studies on soil physics I. The flow of air and water through soils. Journal of Agricultural Research. 4(1), 1–24. [49] Tracy, F. T. 2006. Clean two- and three-dimensional analytical solutions of Richards' equation for testing numerical solvers. Water Resources Research, 42(8), W08503. [50] Brooks, R.H. and Corey, A.T. 1964. Hydraulic Properties of Porous Media and Their Relation to Drainage Design. Transactions of the ASAE, 7(1), 26–28. [51] Lu, N., & Likos, W. J., 2004. Unsaturated soil mechanics. John Wiley and Sons, Inc.

23

Fig. 1 Definition of the elevation head in a slope and illustration of the proposed Chebyshev spectral method.

24

Fig. 2. Test 1: 1D infiltration model of homogeneous unsaturated soil.

25

(a)

(b)

(c)

Fig. 3 Comparison of the calculation accuracy for solving 1D transient flow with different grid sizes and initial conditions using different methods when t = 1 h: (a) h0 = -10 m; (b) h0 = -100,0 m; and (c) h0 = -100,000 m. 26

(a)

(b)

(c)

Fig. 4 Comparison of the computed profiles of the pressure head with different grid sizes at t = 2 h: (a) z  0.5m ; (b) z  0.25m ; and (c) z  0.125m . Test 1 with an initial head of -100,000 m.

27

Fig. 5 Comparison of the calculation results (pressure head) from different methods when h0 = -100,000 m (Test 1).

28

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

29

Fig. 7 Schematic illustration of 2D transient flow.

30

Fig. 8 Boundary condition at the top of the soil.

31

Fig. 9 Comparison of the calculation accuracy for solving 2D transient flow with different grid sizes using different methods when t = 5 h.

32

(a)

(b) Fig. 10 Comparison of results computed using the CSM with the analytical solution for the 2D transient problem at (a) 1 h and (b) 5 h. 33

Fig. 11 Schematic illustration of the homogenous unsaturated slope.

34

(a)

(b) Fig. 12 Fitting curves of the SWCC: (a) sandy soil and (b) silty loam. 35

(a)

(b) Fig. 13 Results comparison with the exact solution for unsaturated soil: (a) sandy soil and (b) silty loam. 36

(a)

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

37

(a)

(b) Fig. 15 Computed results of FS for unsaturated soil: (a) sandy soil and (b) silty loam.

38

Table 1. Parameters of unsaturated soils for analysis. Soil type

K s (m/h)

 g (1/m)

s

r

Sandy soil

6.00 102

1.60 102

0.50

0.11

Silty loam

1.50 103

8.00 103

0.46

0.14

39

Table 2. Variations in the factor of safety (FS) at different depths over the infiltration time for different soils. Conditions

Factor of safety (FS)

Soil type

Depth (m)

0.5 h

1h

1.5 h

2h

Sandy soil

9.263

1.0828

1.0473

1.0315

1.0225

8.247

1.2337

1.1590

1.1240

1.1038

7.27

1.3519

1.2556

1.2071

1.1785

9.263

0.9361

0.8743

0.8423

0.8219

8.247

1.0546

1.0233

0.9928

0.9677

7.27

1.0627

1.0585

1.0480

1.0348

Silty loam

40