Advances in Hydro-Science and Engineering, vol. VIII
Proc. of the 8th International Conference on Hydro-Science and Engineering, Nagoya, Japan, Sept. 9-12, 2008

A NUMERICAL METHOD COMBINING THE IADI METHOD WITH THE
MODIFIED PICARD ITERATION METHOD FOR SOLVING MULTIDIMENSIONAL SATURATED-UNSATURATED FLOW EQUATION
An, H. 1, Ichikawa, Y. 2, Tachikawa, Y. 3 and Shiiba, M.4
1

Student, Department of Urban and Environment Engineering, Kyoto University
Kyotodaigakukatura, Nishikyo-ku, Kyoto, 615-8540, Japan, e-mail: an@hywr.kuciv.kyoto-u.ac.jp
2
Assistant Professor, Department of Urban and Environment Engineering, Kyoto University
Kyotodaigakukatura, Nishikyo-ku, Kyoto, 615-8540, Japan, e-mail: ichikawa@mbox.kudpc.kyoto-u.ac.jp
3
Associate Professor, Department of Urban and Environment Engineering, Kyoto University
Kyotodaigakukatura, Nishikyo-ku, Kyoto, 615-8540, Japan, e-mail: tachikawa@mbox.kudpc.kyoto-u.ac.jp
4
Professor, Department of Urban and Environment Engineering, Kyoto University
Kyotodaigakukatura, Nishikyo-ku, Kyoto, 615-8540, Japan, e-mail: shiiba@mbox.kudpc.kyoto-u.ac.jp
ABSTRACT
Numerical models for simulation of saturated-unsaturated flow are increasingly being
used in water resources evaluation and environmental management. Because of a superiormass balance, the modified Picard iteration method is widely used for solving saturatedunsaturated flow equation. However, the computational cost becomes considerably large
when calculating multi-dimensional flow properties in wide region with finer grid resolution
as it is necessary to make a large matrix calculation to get solution. To reduce thiscomputational cost, in this study, we proposed a numerical method combining the iterative
alternating direction implicit (IADI) method with the modified Picard iteration method forsolving multi-dimensional saturated-unsaturated flow. There is no requirement to make alarge matrix calculation in the process of the proposed method. Two and three dimensionalslope simulations under rainfall conditions are conducted by the proposed method and themodified Picard iteration method. It is found that the proposed method is faster than theconventional modified Picard method while the difference of the simulation results is not
significant. Further it is observed that the new method is more efficient for simulating large
regions with finer resolution.
Keywords: Richards Equation, Modified Picard Method, IADI Method, Saturated-Unsaturated
Flow
1.

INTRODUCTION

Numerical simulation of saturated-unsaturated subsurface flow is being widely used inmany branches of science and engineering including environment engineering, agricultural
engineering, and hydrology. The saturated-unsaturated subsurface flow is governed by theRichards equation (Richards, 1931) which is written as follows∂θ∂K= ∇ ⋅ K (ψ )∇ψ +
∂t
∂z

(1)

Where ψ is the pressure head, θ is the volumetric moisture content, K (ψ ) is hydraulicconductivity, t is time, z is denotes the vertical dimension, assumed positive upward. Also it isassumed that appropriate constitutive relationships between θ and ψ and between K andψ are available. Eq. (1) is called as the mixed form Richards equation including both θ andψ . The mixed form Richards equation has advantages over the ψ -based and θ -based
Richards equation, because the mixed-form is more mass conservative than the ψ -based andθ -based forms (Celia et al. 1990; Clement et al. 1994). Celia et al. (1990) shows that themodified Picard iteration method for the mixed form Richards equation is fully massconserving in the unsaturated zone.
However, the computational cost becomes considerably large when calculating multidimensional flow properties in wide region with finer grid resolution as it is necessary to
make a large matrix calculation to get solution. In this study, we reduce these computationalcosts by combining the iterative alternating direction implicit (IADI) method (Rubin. 1968)
with the modified Picard iteration method. Because the calculation in the process of the IADI
method is conducted implicitly for one direction (dimension) even though the simulation ismulti-dimensional, it can be solved efficiently. To validate the proposed numerical method,
two and three dimensional slope simulations are conducted by the proposed method and themodified Picard iteration method and the result of the two methods are compared.
2.

NUMERICAL METHOD
Modified iteration Picard method
The two dimensional Richards equation is written as follows∂θ ∂ ⎡ ∂ψ ⎤ ∂ ⎡ ⎛ ∂ψ
⎞⎤
= ⎢K
+
K
+
1
⎜
⎟⎥
⎢
∂t ∂x ⎣ ∂x ⎥⎦ ∂z ⎣ ⎝ ∂z⎠⎦(2)where x denotes the horizontal dimension, z denotes the vertical dimension. The modified
Picard iteration method developed by Celia et al. (1990) is based on a fully implicit(backward Euler) time approximation as follows

θ n +1,m +1 − θ nΔt=∂ ⎡ n +1, m ∂ψ n +1,m +1 ⎤ ∂ ⎡ n +1,m ⎛ ∂ψ n +1,m +1 ⎞ ⎤+ ⎢KK+ 1⎟ ⎥⎜∂x ⎢⎣∂x ⎥⎦ ∂z ⎣
∂z
⎝
⎠⎦

(3)

where the superscripts n and m denote time level and iteration level, and the source/sink term
has been ignored for simplicity. The moisture content at the new time step and a new iterationlevel ( θ n +1,m +1 ) is replaced with the truncated Taylor series expansion with respect to ψ ,about the expansion point ψ n +1,m as follows
n +1, mdθθ=θ+(ψ n +1,m +1 −ψ n +1, m ) + 0(δ 2 )(4)dψ
If the higher-order terms are neglected in Eq. (4) and substituting this equation into Eq. (3)gives
n +1, m +1

C n +1,m

ψ n +1, m +1 −ψ n +1,mΔt+n +1, mθ n +1,m − θ nΔt=∂ ⎡ n +1,m ∂ψ n +1,m +1 ⎤ ∂ ⎡ n +1, m ⎛ ∂ψ n +1,m +1 ⎞ ⎤+ ⎢KK+ 1⎟ ⎥ (5)⎜∂x ⎣⎢∂x ⎥⎦ ∂z ⎣∂z⎝⎠⎦where C = dθ / dψ is the specific moisture capacity function. Eq. (5) can be rewritten in theterms of the increment in iteration δ m = ψ n +1,m +1 −ψ n +1, m as follows
C n +1, m m ∂ ⎡ n +1,m ∂δ m ⎤ ∂ ⎡ n +1,m ∂δ m ⎤δ + ⎢K+KΔt∂x ⎣∂x ⎥⎦ ∂z ⎢⎣
∂z ⎥⎦
=−

θ n +1,m − θ nΔt∂ ⎡∂ψ n +1, m ⎤ ∂ ⎡ n +1,m ⎛ ∂ψ n +1, m ⎞ ⎤+ ⎢ K n +1, m+ ⎢K+ 1⎟ ⎥⎜∂x ⎣∂x ⎥⎦ ∂z ⎣
∂
z
⎝
⎠⎦

(6)

Eq. (6) defines the modified Picard iteration method of Celia et al. (1990). δ m is theunknown dependent variable and this simultaneous equation is able to be solved by making amatrix calculation. In case of multi-dimensional simulation for large region with fine grid
resolution, a large matrix calculation is required in the process of calculating simultaneous
equation.
The iterative process of Eq. (6) continues until the difference between the calculated
values of the moisture content of two successive iteration levels becomes less than thetolerance δθ (Huang et al. 1996) as follows

θ n +1,m +1 − θ n +1, m ≤ δθ(7)

In this study, δθ = 0.00001 is used. Huang et al. (1996) showed that this convergencecriterion based on moisture content is more efficient in a computation than based on pressure
head especially when the soil hydraulic characteristics are highly nonlinear.
IADI Method and Its Coupling with the Modified Picard Method
Iterative alternating direction implicit (IADI) method of Rubin (1968) is an algorithmfor ψ -based form Richards equation as follows
C (ψ )∂ψ∂ ⎡ ∂ψ ⎤ ∂ ⎡ ⎛ ∂ψ⎞⎤
= ⎢K
+ ⎢K ⎜+ 1⎟ ⎥⎥∂t ∂x ⎣ ∂x ⎦ ∂z ⎣ ⎝ ∂z⎠⎦

(8)

The linearized expressions of Eq. (8) based on the IADI method is written as follows
Cin, +j 1,2 m

ψ in, +j 1,2 m +1 −ψ in, jΔt+ I m K n (ψ in, +j 1,2 m +1 −ψ in, +j 1,2 m )⎞⎤∂ ⎡ ∂ψ n +1,2 m +1 ⎤∂ ⎡ n ⎛ ∂ψ n +1,2 m= ⎢K n++ 1⎟ ⎥
⎢K ⎜
⎥
∂x ⎣
∂x
⎦ i , j ∂z ⎣ ⎝ ∂z⎠⎦i, jn +1,2 mi, jCψ in, +j 1,2 m + 2 −ψ in, jΔt(9)+ I m K n (ψ in, +j 1,2 m + 2 −ψ in, +j 1,2 m +1 )∂ ⎡ ∂ψ n +1,2 m +1 ⎤∂ ⎡ n ⎛ ∂ψ n +1,2 m + 2 ⎞ ⎤= ⎢K n++ 1⎟ ⎥
⎢K ⎜
⎥
∂x ⎣
∂x
∂z
⎠⎦i, j⎦ i , j ∂z ⎣ ⎝(10)

K in, j = K in−1/ 2, j + K in+1/ 2, j + K in, j −1/ 2 + K in, j +1/ 2(11)

where subscript i and j means the spatial coordinates of the point in x and z axis respectively,
and I m is the iteration parameter, computed for each iteration. In this study, I m = 0.55m isused according to Weeks (2004). Because Eqs (9) and (10) are one-directional implicit forms,
there is no requirement to make a large matrix calculation for solving the equations eventhough in a multi-dimension.
The Richards equation can be efficiently solved by using IADI method. But as thisscheme is an algorithm for the ψ -based form Richards equation, it is often faces to a
problem of mass conservation error. To take advantage of two numerical methods, IADI
method and the modified Picard iteration method, we propose a numerical method combining
both two methods. The proposed method applied to Eq. (2) is written as follows

ψ in, +j 1,2 m +1 −ψ in, +j 1,2 m

θin, +j 1,2 m − θ in, j+ I m K n (ψ in, +j 1,2 m +1 −ψ in, +j 1,2 m )
Δt
⎞⎤
∂ ⎡
∂ψ n +1,2 m +1 ⎤∂ ⎡ n +1,2 m ⎛ ∂ψ n +1,2 m= ⎢ K n +1,2 m++ 1⎟ ⎥
⎢K
⎜
⎥
∂x ⎣
∂x
⎦ i , j ∂z ⎣⎝ ∂z⎠⎦i, jn +1,2 m
i, j

C

Δt

Cin, +j 1,2 m +1+

ψ in, +j 1,2 m + 2 −ψ in, +j 1,2 m +1 θ in, +j 1,2 m +1 − θ in, j
+

Δt

∂ ⎡
∂ψ
= ⎢ K n +1,2 m∂x ⎣∂x

n +1,2 m +1

Δt

(12)

+ I m K n (ψ in, +j 1,2 m + 2 −ψ in, +j 1,2 m +1 )⎤∂ ⎡ n +1,2 m ⎛ ∂ψ n +1,2 m + 2 ⎞ ⎤++ 1⎟ ⎥⎢K⎜⎥∂z⎦ i , j ∂z ⎣⎝⎠⎦i, j(13)

These equations are rewritten as follows
n +1,2 m
i, j

C

=−

δ 2m

θin, +j 1,2 m − θin, jΔtCin, +j 1,2 m +1=−

∂ ⎡ n +1,2 m ∂δ 2 m ⎤+ ⎢K+ I m K nδ 2 m⎥
Δt ∂x ⎣∂x ⎦ i , j

δ 2 m +1
Δt

+

⎞⎤
∂ ⎡
∂ψ n +1,2 m ⎤∂ ⎡ n +1,2 m ⎛ ∂ψ n +1,2 m+ ⎢ K n +1,2 m++ 1⎟ ⎥⎢K⎜⎥∂x ⎣∂x ⎦ i , j ∂z ⎣⎠⎦i, j⎝ ∂z∂ ⎡ n +1,2 m ∂δ 2 m +1 ⎤K
+ I m K nδ 2 m +1⎢⎥
∂z ⎣∂z ⎦ i , jθin, +j 1,2 m +1 − θ in, jΔt∂ ⎡∂ψ n +1,2 m +1 ⎤∂ ⎡ n +1,2 m ⎛ ∂ψ n +1,2 m +1 ⎞ ⎤+ ⎢ K n +1,2 m++ 1⎟ ⎥
⎢K
⎜
⎥
∂x ⎣
∂x
∂z
⎦ i , j ∂z ⎣⎝⎠⎦i, j(14)(15)

Eqs (14) and (15) are one-directional implicit forms and based on mixed-form Richards
equation. δ 2m and δ 2 m +1 are the unknown dependent variables and the convergencecriterion is as follows

θ n +1,2 m + 2 − θ n +1,2 m ≤ δθ(16)

The numerical method using Eqs (14) and (15) was confirmed by conducting two dimensionalslope simulation.

Modifying IADI method
Extending Eqs (14) and (15) for a three-dimensional field, it is written as followsmCin, +j ,1,3
k

=−

δ 3m

+

Δt

∂ ⎡ n +1,3m ∂δ 3m ⎤K
+ I m K nδ 3 m∂x ⎢⎣∂x ⎥⎦ i , j ,km− θin, j ,kθin, +j ,1,3
k

Δt

+

⎞⎤
∂ ⎡ n +1,3m ∂ψ n +1,3m ⎤∂ ⎡ n +1,3m ∂ψ n +1,3m ⎤∂ ⎡ n +1,3m ⎛ ∂ψ n +1,3m+++ 1⎟ ⎥
K
K
⎢K
⎜
⎢
⎥
⎢
⎥
∂x ⎣
∂x ⎦ i , j ,k ∂y ⎣∂y ⎦ i , j ,k ∂z ⎣⎝ ∂z
⎠ ⎦ i , j ,k

(17)
m +1
Cin, +j ,1,3
k

=−

δ

3 m +1Δt+⎤∂ ⎡ n +1,3m ∂δ+ I m K nδ 3m+1
K⎢∂y ⎣∂y ⎥⎦ i , j , k
3 m +1m +1− θin, j , kθin, +j ,1,3
k

Δt

+

∂ ⎡ n +1,3m+1 ∂ψ n +1,3m +1 ⎤∂ ⎡ n +1,3m +1 ∂ψ n +1,3m +1 ⎤∂ ⎡ n +1,3m +1 ⎛ ∂ψ n +1,3m +1 ⎞ ⎤+++ 1⎟ ⎥
K
K
⎢K
⎜
⎥
⎢
⎥
∂x ⎢⎣
∂x
∂y
∂z
⎦ i , j ,k ∂y ⎣⎦ i , j ,k ∂z ⎣⎝⎠ ⎦ i , j ,k

(18)
m+2
Cin, +j ,1,3
k

=−

δ

3m + 2Δt+∂ ⎡ n +1,3m+ 2 ⎛ ∂δ
⎢K
⎜
∂z ⎣
⎝ ∂z

m+ 2
− θ in, j ,kθin, +j ,1,3kΔt

3m + 2+⎞⎤+ 1⎟ ⎥ + I m K nδ 3m+ 2⎠ ⎦ i , j ,k⎛ ∂ψ n +1,3m + 2 ⎞ ⎤∂ ⎡ n +1,3m + 2 ∂ψ n +1,3m + 2 ⎤∂ ⎡∂ψ n +1,3m + 2 ⎤∂ ⎡+ ⎢ K n +1,3m + 2+ ⎢ K n +1,3m + 2 ⎜+ 1⎟ ⎥
K
⎢
⎥
⎥
∂x ⎣
∂x
∂y
∂z
⎦ i , j ,k ∂y ⎣⎦ i , j ,k ∂z ⎣⎝⎠ ⎦ i , j ,k(19)
where subscript i, j, and k are spatial coordinates of the point in x, y, and z axis respectively.
We conducted three dimensional slope simulation using Eqs (17)~(19). However, thecalculation using these equations was unstable and the iterative calculation did not converge.
The reason of this instability has not been made clear, but we supposed that the terms
calculated implicitly are not enough to solve the equation in a three dimension simulation.
The two and three dimensional equations, Eqs (14)~(15), (17)~(19), are one-directional
implicit forms. The half of the direction are implicitly calculated in the two-dimensional
equations, Eqs (14) and (15). On the other hand, the one of three of the directions are
implicitly calculated in the three-dimensional equations, Eqs (17)~(19). Therefore, it is
supposed to be needed to calculate implicitly over the half of the direction terms forstabilizing the iterative process.
To stabilize the iterative calculation in a three-dimensional simulation, we modified
the IADI method part of the proposed method. It is assumed that Eq. (17) is calculated
according to the order shown in Fig.1 . Then, when the red region of Fig. 1 is being calculated,the values of the next iteration step in the blue region were already solved. Hence, we can use
these values for calculating the red region. According to this idea, we modified Eq. (17) as

Figure 1 the order of the calculation for equation (17) in modified IADI method

n +1,3 m
i , j ,k

C

δ 3m
Δt

+

m
θin, +j ,1,3
− θin, j ,k
k

Δt

1
m +1
m +1
⎡ K in−+1/1,32,mj , kψ in−+1,1,3j , km +1 + K in−+1/1,32,mj ,kψ in, +j ,1,3
+ K in++1/1,32,mj ,kψ in, +j ,1,3
+ K in++1/1,32,mj ,kψ in++1,1,3j , km+1 ⎤⎦kk
Δx 2 ⎣
1m +1
m+ 2 ⎡⎣ K in, +j −1,31/m2,kψ in, +j −1,31,km +1 + K in, +j −1,31/m2, kψ in, +j ,1,3
+ K in, +j +1,31/m2, kψ in, +j ,1,3
+ K in, +j +1,31/m2,kψ in, +j −1,31,km ⎤⎦
k
k
Δy
1
m
n +1,3 mm
n +1,3 m +1
m
n +1,3 m +1
m
n +1,3 m⎤+ Kin, +j ,1,3
+ 2 ⎡⎣ K in, +j ,1,3
+ K in, +j ,1,3
+ K in, +j ,1,3
k +1/ 2ψ i , j , k +1 ⎦k −1/ 2ψ i , j , k −1
k −1/ 2ψ i , j , k
k +1/ 2ψ i , j , kΔz1mn +1,3 m⎡ K in, +j ,1,3⎤+
k +1/ 2 − K i , j , k −1/ 2 ⎦
⎣
Δz
=

(20)
This equations is rewritten as follows

δ i3,mj ,k1⎡ K n +1,3m δ 3m + K in−+1/1,32,mj ,k δ i3, mj ,k + K in++1/1,32,mj ,k δ i3, mj , k + K in++1/1,32,mj , k δ i3+m1, j ,k ⎤⎦2 ⎣ i −1/ 2, j , k i −1, j , kΔtΔxn +1,3 m
3m
m
3m
K i , j −1/ 2,k δ i , j ,k K in, +j ,1,3
k −1/ 2δ i , j , k
+
+
Δy 2
Δz 2
n +1,3 m
i , j ,k

C

=−

+

m
− θin, j ,kθin, +j ,1,3
k

Δt

1
m
m
⎡ K in−+1/1,32,mj ,kψ in−+1,1,3j , km + K in−+1/1,32,mj , kψ in, +j ,1,3
+ K in++1/1,32,mj ,kψ in, +j ,1,3
+ K in++1/1,32,mj ,kψ in++1,1,3j ,km ⎤⎦kk
Δx 2 ⎣1m
m+ K in, +j +1,31/m2,kψ in, +j ,1,3
+ K in, +j +1,31/m2, kψ in, +j −1,31,km ⎤⎦+ 2 ⎡⎣ K in, +j −1,31/m2,kψ in, +j −1,31,km +1 + K in, +j −1,31/m2,kψ in, +j ,1,3
k
k
Δy
1
m
n +1,3 m +1
m
n +1,3 mm
n +1,3 mm
n +1,3 m⎤+ 2 ⎡⎣ K in, +j ,1,3
+ K in, +j ,1,3
+ K in, +j ,1,3
+ K in, +j ,1,3
k −1/ 2ψ i , j , k −1
k −1/ 2ψ i , j , k
k +1/ 2ψ i , j , k
k +1/ 2ψ i , j , k +1 ⎦Δz1m
n +1,3 m⎡⎣ K in, +j ,1,3⎤+
k +1/ 2 − K i , j , k −1/ 2 ⎦
Δz
+

(21)
Eq. (20) is implicit for not only x direction but also the half of y and z directions. The terms of
m +1
ψ in, +j −1,31,km+1 and ψ in, +j ,1,3
in the right side of Eq. (21) are already calculated values and δ 3m is
k −1
the only unknown dependent variable in Eq. (21). Therefore, the equation implicit for x
direction and the half of y and z directions is solved in the process of this method but thisequation can be solved by calculating the equation implicit for only one direction. The
equation of y and z directions, Eqs (18) and (19), are also modified in the same way as theequation of x direction.
This numerical method combining the modified IADI method with modified Picard
iteration method is a newly developed method by this study.
3.

NUMERICAL EXPERIMENTS
Two-dimensional simulations

To validate the proposed computation method, the two-dimensional slope simulations
are conducted by using the modified Picard iteration method and the proposed method. The

simulation condition is shown in Fig. 2. Soil slope, 10m long, 1m thickness and 20 degree
inclined, is concerned. A constant water level is maintained at the lower end and the bottomside and the upper end are no flux boundaries. The slope surface is under rainfall condition
and Fig. 2(b) shows the rainfall intensity. Table 1 shows the four grid size conditions used in
the simulation. The time step is 10 sec and the simulation length is 10 hours. The soil water
retention and hydraulic conductivity were described as follows⎛ ψ′⎞θ − θr ⎛ ψ ′ ⎞=⎜+ 1⎟ exp ⎜ − ⎟ ,θs − θr ⎝ψ 0 ⎠⎝ ψ0 ⎠⎛ θ − θr ⎞K = Ks ⎜⎟⎝ θs − θr ⎠

⎧ψ
ψ′= ⎨

⎩0

(ψ < 0)(ψ ≥ 0)

(22)

β

(23)

where θ s and θ r are the saturated and residual water contents respectively; K s is thesaturated hydraulic conductivity; ψ 0 is air entry pore-water pressure head; β is theparameter determining the shape of unsaturated hydraulic conductivity curve. In this study,θ s =0.5, θ r =0.1, ψ 0 =-0.25m, K s = 0.00001m/s, β =3.0 are used.

Figure 2 Simulation condition: (a) slope area, (b) rainfall intensity
Table 1 Grid size conditions for two-dimensional simulation

Case1
Case2
Case3
Case4

+x
0.1 m
0.05 m0.05 m0.25 m+z
0.05 m0.05 m0.25 m0.25 mNumber of mesh
25000
50000
100000200000Figure 3 shows the error which is defined as follows
Error =

θ MP − θCo
θ MP

(24)

where θ MP is the moisture content calculated by the modified Picard iteration method, θ Co
is the moisture content calculated by the proposed method. Figure 3(a) is the average error
and Figure 3(b) is the maximum error at all grids points. The peak value of the average errors

are around 5.0 × 10−5 and the peak value of the maximum errors are around 1.4 × 10−3 . Bothvalues are small enough to be neglected.
Figure 4(a) shows computation times of the modified Picard iteration method and of
the proposed method for each case. Figure 4(b) shows the ratio of the computation time of the
proposed method to that of the modified Picard iteration method. In this simulation, theproposed method is faster from 1.5 to 5 times than the modified Picard iteration method.

Figure 3 Error of moisture content: (a) average value, (b) maximum value in all grids

Figure 4 Computation times in two-dimensional simulations
Three-dimensional simulations
The three-dimensional simulation condition is shown in Fig. 5. Soil slope, 10m long,
5m wide, 1m thickness and 20 degree inclined, is concerned. The half of the surface is underrainfall condition so that the subsurface flow becomes three-dimensional flow. The Rainfall
intensity is shown in Fig. 2(b). A constant water level is maintained at the lower end and theboth lateral sides, the bottom side, and the upper end are no flux boundaries. Table 3 shows
the four grid size conditions used in the simulation. The time step is 10 sec and the simulation
length is 10 hours. The soil water retention and hydraulic conductivity were same in the twodimensional simulations.
Table 3 Grid size conditions for three-dimensional simulation

Case1
Case2
Case3
Case4

+x
0.2 m
0.2 m
0.1 m
0.1 m

+y
0.2 m
0.1 m
0.1 m
0.1 m

+z
0.05 m0.05 m0.05 m
0.025 mNumber of mesh
25000
50000
100000200000Figure 5 Simulation condition
Figure 6 shows the error which is defined as Eq. (24). The peak value of the average
errors are around 2.5 × 10−4 and the peak value of the maximum errors are around 6.5 × 10−3 .
Both values are small enough to be neglected.
Figure 7(a) shows computation times of the modified Picard iteration method and of
the proposed method for each case. Figure 7(b) shows the ratio of the computation time of the
combining method to that of the modified Picard iteration method. In this simulation, theproposed method is faster from 6 to 8 times than the modified Picard iteration method.

Figure 6 Error of moisture content: (a) average value, (b) maximum value in all grids

Figure 7 Computation times in three-dimensional simulations
4.

CONCLUSION

This study proposed a new computation method combining the iterative alternating
direction implicit (IADI) method with the modified Picard iteration method for solving multidimensional saturated-unsaturated flow. There is no requirement to make a large matrix
calculation in the process of the proposed method.
To validate the proposed method, the two and three dimensional slope simulations
were conducted. The proposed method was stable in the two-dimensional simulations. But inthe three-dimensional simulation, it was unstable and iterative calculation did not converge.
To stabilize the iterative calculation in the three-dimensional simulations, we modified the
IADI method part of the proposed method. This method combining the modified IADI

method with the modified Picard iteration method is applicable in the two and three
dimensional simulations.
The two and three dimensional slope simulations under rainfall conditions wereconducted by the proposed method and the modified Picard iteration method. It is found thatthe proposed method is faster than the conventional modified Picard method while thedifference of the simulation results is not significant. Further it is observed that the new
method is more efficient for simulating with finer resolution.
REFERENCES
Celia, M.L., Bouloutas, E.T., Zarba R.L. (1990), A General Mass-Conservative Numerical
Solution for the Unsaturated flow Equation, Water Resources Research, 26(7), pp.14831496.
Clement, T.P., Wise, W.R., Molz, F.J. (1994), A Physically Based, Two-dimensional, FiniteDifference Algorithm for Modeling Variably saturated Flow, Journal of Hydrology, 161,
pp.71-90.
Huang, K., Mohanty, B.P., van Genuchten M.Th. (1996), A New Convergence Criterion for
the Modified Picard Iteration Method to Solve the Variably Saturated Flow Equation,
Journal of Hydrology, 178, pp.69-91.
Richards, L.A. (1931), Capillary conduction of liquids in porous mediums, Physics, 1,
pp.318-333.
Rubin, J. (1968), Theoretical Analysis of Two-Dimensional, Transient Flow of Water in
Unsaturated and Partly Unsaturated Soils, Soil Science Society of America, 32(5), pp.607615.
Weeks, S.W., Sander, G.C., Braddock, R.D., Matthews, C.J. (2004), Saturated andUnsaturated water flow in inclined porous media, Environmental Modeling andAssessment, 9, pp.91-102.