A modified Picard iteration scheme for overcoming numerical difficulties of simulating infiltration into dry soil

Journal of Hydrology 551 (2017) 56–69

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Research papers

A modified Picard iteration scheme for overcoming numerical difficulties of simulating infiltration into dry soil Yuanyuan Zha a,b, Jinzhong Yang a, Lihe Yin b, Yonggen Zhang c, Wenzhi Zeng a,d, Liangsheng Shi a,⇑ a

State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan, Hubei 430072, China Key Laboratory for Groundwater and Ecology in Arid and Semi-Arid Areas, China Geological Survey, Xi’an, Shaanxi 710054, China c Department of Hydrology and Atmospheric Sciences, University of Arizona, Tucson, AZ 85721, USA d Crop Science Group, Institute of Crop Science and Resource Conservation (INRES), University of Bonn, Katzenburgweg 5, D-53115 Bonn, Germany b

a r t i c l e

i n f o

Article history: Received 28 January 2017 Received in revised form 24 May 2017 Accepted 25 May 2017 Available online 27 May 2017 This manuscript was handled by Corrado Corradini, Editor-in-Chief, with the assistance of Renato Morbidelli, Associate Editor Keywords: Variably saturated flow Richards’ equation Model robustness Picard iteration Infiltration

a b s t r a c t Numerical models based on Richards’ equation are often employed to simulate the soil water dynamics. Among them, those Picard iteration models which use the head as primary variable are widely adopted due to their simplicity and capability for handling partially saturated flow conditions. However, it is wellknown that those models are prone to convergence failure in some unfavorable flow conditions, especially when simulating infiltration into initially dry soils. Here we analyze the reasons that give rise to the numerical difficulty. Moreover, several modifications to the mass-conservative Picard iteration method are proposed so that numerical difficulty is avoided in these unfavorable flow conditions. Our proposed modifications do not degrade the simulated results, while they lead to more robust convergence performances and cost-effective simulations. Ó 2017 Elsevier B.V. All rights reserved.

1. Introduction Modeling variably saturated flow in the vadose zone is commonly based on Richards’ equation, which is derived based on the combination of flow continuity and Darcy’s law. With specific descriptions of the soil hydraulic functions, analytical solutions of Richards’ equation can be obtained for those problems with simple geometric condition as well as initial and boundary conditions (Menziani et al., 2007; Srivastava and Yeh, 1991; Tracy, 2006). To deal with complex problems, numerical models are more general. To solve Richards’ equation numerically, the soil water retention curve is introduced to eliminate one of the two dependent variables, i.e., pressure head and soil moisture. Accordingly, depending on which one is used as the primary variable, numerical algorithms can be classified into pressure head-based form, moisture-based form, or a combination of these two (i.e., the socalled primary switching technique). The primary switching technique (Diersch and Perrochet, 1999; Forsyth et al., 1995; Sadegh Zadeh, 2011), which uses soil moisture as primary variable for

⇑ Corresponding author. E-mail addresses: [email protected], [email protected] (L. Shi). http://dx.doi.org/10.1016/j.jhydrol.2017.05.053 0022-1694/Ó 2017 Elsevier B.V. All rights reserved.

unsaturated nodes but pressure head for saturated nodes, was found to yield rapid convergence when the soil is dry (Ross, 2003; Hassane Maina and Ackerer, 2016). However, the nonsmooth transition between alternative primary variables could lead to unrealistic results (e.g., moisture content is greater than 1) and this technique relies on Newton solver, which is quite demanding from an implementation point of review (Krabbenhøft, 2007). In contrast, the algorithms that only use one primary variable (either pressure head or soil moisture) is free from the nonsmooth transition issue, and it can be implemented with Picard iteration solver, which is more straightforward since it does not require the calculation of Jacobian matrix. However, due to the nonlinear relationship between head and soil moisture, the transformed formulas (i.e., calculating Darcian flux in terms of soil moisture or measuring mass change using head) have prominent flaws. For instance, although the moisture-based form equation is inherently mass conservative, the hydraulic gradient cannot be correctly represented in the form of soil moisture either in saturated or in heterogeneous soils (Crevoisier et al., 2009). Despite there are several treatments (e.g., Hills et al., 1989; Kirkland et al., 1992; Matthews et al., 2004; Zha et al., 2013b) to fix this problem in the heterogeneous soil, this form is inappropriate to simulate unsaturated-saturated flow and thus has limited its applications.

57

Y. Zha et al. / Journal of Hydrology 551 (2017) 56–69

In contrast, the head-form algorithm is often employed for variably saturated flow conditions and flow in heterogeneous soils. The Picard iteration solvers using the head as primary variable are widely used since a number of popular codes/software packages (e.g., Šimu˚nek et al., 2009; van Dam and Feddes, 2000) are available. However, this algorithm also has some drawbacks that hinder its popularity. First, temporal discretization of the storage term can introduce significant mass balance errors (Hao et al., 2005; Lai and Ogden, 2015). Celia et al. (1990) presented a numerical algorithm where the soil moisture at current iteration level, obtained by the newest pressure head, is used to measure the mass change. In essence, Celia method utilizes the truncated Taylor series of soil moisture with respect to pressure head about the iteration point. Since it can maintain the global mass balance as long as the solution converges with iterations (Kosugi, 2008), Celia scheme implemented in Picard iteration model has been widely adopted although there are other approaches (e.g., Rathfelder et al., 1994). Later on, to further enhance the mass conservation performance, it is suggested to apply Celia algorithm only if the change of pressure head between time steps is larger than the threshold value (e.g., 3–5 cm) (Hao et al., 2005) or the nodal head value is less than the threshold value (e.g.,
2.5 cm) (Sadegh Zadeh, 2011). Nevertheless, the mass balance error of Celia algorithm or its variants relies on the converging error (the head difference between adjacent iterations). Second, when simulating flow into initially dry soil, the head gradient at the wetting front is extremely large, leading to large computational efforts, numerical instability, and inaccuracy. The error of the moisture profile can be effectively reduced by the introduction of adaptive time stepping (Kavetski et al., 2001) and mesh refinement (Miller et al., 2006) with additional computational cost. On the other hand, the numerical divergence problem is more serious since it leads to the termination of simulation. Huang et al. (1996) proposed a new convergence criterion for the Celia et al. (1990) algorithm. In their approach, convergence rests on the difference of soil moisture or pressure head at two adjacent iteration levels when the node is unsaturated and saturated, respectively. The numerical experiments showed that the new criterion can both reduce the computational costs and the chance of divergence. It should be noted that their simulations were restricted to the constant-head upper boundary. On the contrary, simulation interrupt due to divergence is often confronted when the soil is subject to alternate wetting-drying flux conditions (i.e., precipitation and evaporation). Although the inefficiency of headform algorithm for infiltration into dry soil has been recognized (e.g., Forsyth et al., 1995; Crevoisier et al., 2009; Varado et al., 2006a,b, Zha et al., 2016, 2013a), a thorough investigation of the divergence phenomenon is still lacking. The origin of the difficulty, the impact of different soil hydraulic functions as well as flow conditions for this problem require comprehensive assessment. Even though the Picard iteration solver with Celia et al. (1990) algorithm can achieve good mass balance and accurate moisture profile with affordable cost, these good performances can be destroyed when the flow has difficulty in convergence under these unfavorable flow conditions. Development of a more stable numerical scheme is urgent for researchers who focus on variably saturated flow modeling. Thus, to further promote the use of Picard iteration solvers in widely-used codes/software packages (e.g., Šimu˚nek et al., 2009; van Dam and Feddes, 2000), we extensively assess the robustness of the mass-conservation Picard iteration scheme proposed by Celia et al. (1990). The reasons for the numerical instability are thoroughly analyzed for the scenarios with different combinations of initial conditions, infiltration intensity, mesh size, and soil hydraulic properties (Section 2). Starting from the Celia

scheme, we propose to modify the first-cut head value in the iteration procedure, so that the requirement of time step size originating from the nonlinearity can be relaxed (Section 3). A remarkable increase in model robustness using the modified algorithm is attained, as demonstrated by several numerical experiments (Section 4). Summary and conclusions are presented in Section 5. 2. Numerical formulation 2.1. Governing equations For three-dimensional flow in unsaturated-saturated soils, Richards’ equation is the combination of the mass conservation equation and the Darcy’s Law:

@h @h þ bls ¼
r  q @t @t

ð1Þ

q ¼
KðhÞrH

ð2Þ 3
3

where t is time [T]; h is volumetric moisture content [L L ]; H [L] is the hydraulic head equal to pressure head h plus elevation head z; q is Darcian flux vector [LT
1] and K is the unsaturated hydraulic conductivity [LT
1]. The saturation index b is equal to one if the medium is saturated and zero if the soil is unsaturated (Mao et al., 2011). ls [L
1] is the specific storage. When only one-dimensional vertical flux is concerned, q and r  q in Eq. (2) are simplified as qz =
K(h)(1 + oh/oz) and oqz/oz. Substituting Eq. (2) into Eq. (1) leads to the mixed form Richards’ equation. Transformation of the storage term into head leads to head-based form Richards’ equation, while calculating Darcian flux from soil moisture results in moisture-based form Richards’ equation. The soil water retention curve as well as the unsaturated hydraulic conductivity must be provided to describe the relationships between the unknown soil hydraulic parameters and the unknown primary variable. Moreover, initial and boundary conditions are required to obtain the solution. 2.2. Numerical scheme Many approaches can be employed to discretize the spatial derivatives in Richards’ equation, such as finite element method (Šimu˚nek et al., 2009), finite difference method (Dogan et al., 2005), finite volume method (Caviedes-Voullième et al., 2013), or finite analytical method (Zhang et al., 2015). Since it is well documented and it is beyond the scope of this study, the spatial discretization formulation is not discussed here and we only focus on the temporal discretization scheme. As discussed in Kavetski (2002), an efficient time-stepping scheme as well as a reasonable temporal derivative discretization are also important factors of a cost-effective model. The original Picard iteration discretizes əh/ət as follows: jþ1;kþ1

@h @h h ¼C  C jþ1;k @t @t

Dt


h

j

ð3Þ

where Dt is the time step [T], C is the soil water capacity [L
1], and the superscript j and k represent the time level and iteration level, respectively. Although mathematically əh/ət = Cəh/ət, numerically Cj+1,k(hj+1,k+1
hj) is not equal to (hj+1,k+1
hj) unless the change of head (hj+1,k+1
hj) is infinitesimal (Celia et al., 1990). Thus, this discretization could lead to large mass balance error even if the model converges (|hj+1,k+1
hj+1,k|  0). Celia et al. (1990) proposed a modified Picard iteration algorithm (hereinafter referred to as CI algorithm) to ensure the mass balance during the temporal discretization,

58

Y. Zha et al. / Journal of Hydrology 551 (2017) 56–69

convergence criteria for unsaturated (h < 0) and saturated (h  0) nodes are,

@h hjþ1;kþ1
hjþ1;k hjþ1;k
h j ¼ þ @t Dt Dt jþ1;kþ1

¼ C jþ1;k

h


h Dt

jþ1;k

þ

hjþ1;k
h j Dt

ð4Þ

When it converges, the first term on the right side of Eq. (4) will approach 0. At the first iteration (k = 0), the solution from the previous time step (or initial condition) is usually used as the initial guess, i.e., hj ? hj+1,0. Once a solution of hj+1,k is obtained from solving the equation set, the corresponding soil moisture hj+1,k can be calculated immediately via the soil water retention curve. After that, the results hj+1,k and hj+1,k are substituted into Eq. (4) for solving the head at new iteration level (hj+1,k+1). The iteration level k marches (k k + 1) until the solution converges, i.e., the predefined convergence criterion is satisfied (Huang et al., 1996). The CI algorithm has been widely used due to its simplicity (e.g., software HYDRUS (Šimu˚nek et al., 2009) and code SWAP (van Dam and Feddes, 2000)).

max jhjþ1;kþ1
hjþ1;k j < eh ; i for all unsaturated nodes i i jþ1;kþ1

max jhi

jþ1;k


hi

ð5Þ

j < eh ; i for all saturated nodes

ð6Þ

Unless otherwise noted, eh is given as 10 cm cm and eh is given as 0.1 cm for the numerical examples in Section 4. The current time-step size is adjusted according to the convergence performance at previous time step. In order to reduce the chance of convergence failure, the program recovers the head values before the iteration loop (hj ? hj+1,0) and a recalculation at current time level with a much smaller Dt is conducted if the maximum number of iterations is reached during one time step, as done in some popular software, e.g., HYDRUS-1D (Šimu˚nek et al., 2009) and VSAFT3 (Srivastava and Yeh, 1992). This procedure is recursive until Dt is smaller than a user-specified minimum time step Dtmin.
4

3


3

2.4. Analysis of the numerical difficulties 2.3. Convergence criteria and time stepping scheme Dynamic time stepping can be easily handled for the iterative Newton or Picard schemes (D’Haese et al., 2007). The time step size can be increased if convergence at the current time level is achieved in very few iterations (e.g., iteration number nit < 3), and decreased if the solution converges slowly (e.g., nit > 7). The

2.4.1. Reason for numerical difficulties It is well known that CI method still suffers from divergence for problems involving infiltration into dry soils. Fig. 1 presents soil water capacity and hydraulic conductivity curves for loam soil (the van Genuchten parameters (van Genuchten, 1980) from Carsel and Parrish (1988) are listed in Table 1). The dots on the

Fig. 1. Change of intermediate capacity C (a) and unsaturated hydraulic conductivity K (b) of surface node during iterations for simulating infiltration into initially dry soil. The index k represents iteration level and hj+1,* is the convergent solution at time step j + 1.

59

Y. Zha et al. / Journal of Hydrology 551 (2017) 56–69 Table 1 The soil hydraulic parameters from Carsel and Parrish (1988) and related hwet and hdry values used in the numerical cases. Soil type

hr (cm
3 cm
3)

hs (cm
3 cm
3)

a (cm
1)

n

Ks (cm d
1)

hwet (cm)

hdry (cm)

Sand Loamy sand Sandy loam Loam Silt Silt loam Sandy clay loam Clay loam Silty clay loam Sandy clay Silty clay Clay

0.045 0.057 0.065 0.078 0.034 0.067 0.1 0.095 0.089 0.1 0.07 0.068

0.43 0.41 0.41 0.43 0.46 0.45 0.39 0.41 0.43 0.38 0.36 0.38

0.145 0.124 0.075 0.036 0.016 0.02 0.059 0.019 0.01 0.027 0.005 0.008

2.68 2.28 1.89 1.56 1.37 1.41 1.48 1.31 1.23 1.23 1.09 1.09

712.8 350.2 106.1 24.96 6 10.38 31.44 6.24 1.68 2.88 0.48 4.8


2.6
2.3
2.3
2.1
1.7
1.8
0.9
0.8
0.5
0.2 0.0 0.0


209
355
801
1865
3527
2942
1577
2823
3426
2441
1949
2075

curves indicate the coordinates (h, C) and (h, K) during different iteration levels (denoted by superscript k). According to CI scheme, the discrete water balance equation for the surface node is (assume 1-D flow), jþ1;kþ1

C jþ1;k 1

h1

jþ1;k


h1 Dt

jþ1;kþ1

hjþ1;k
h1j h 1
bjþ1;k ls 1 1 Dt 1 jþ1;kþ1 ðI
q1þ1=2 Þ þ Dz1

¼


Dt

j


h1 ð7Þ

where I is the infiltration rate, Dz1 = 0.5 cm is the controlled thickjþ1;kþ1

jþ1;k ness of the surface node 1, qjþ1;kþ1 1þ1=2 ¼ K 1þ1=2 ðh1

ere, since the boundary conditions constrain the head change. However, the flux boundary is commonly used since the precipitation intensity can be easily obtained rather than the pressure head measurement.

jþ1;kþ1


h2

Þ=Dz1þ1=2

is the Darcian flux (positive downward) between node 1 and node 2, and Dz1+1/2 = 1 cm is the distance between the two nodes. Initially, the surface node is relatively dry (hj 
3000 cm and hj  0.1), while both Cj (10
5 cm
1 in Fig. 1) and Kj (10
4 cm/d in Fig. 1) are very small. At the first iteration (k = 0), the calculated Darcian flux qj+1,1 should also be very small since it is determined by Kj+1,0 = Kj (lagged by one iteration level as done in Picard iteration). Furthermore, hj+1,0 = hj and hj+1,0 = hj are used as initial guesses. Under this situation, the change of head during one iteration (i.e., hj+1,0+1–hj+1,0) is tremendous and the obtained head hj+1,0+1 often approaches or even exceed 0 (Fig. 1) even though Dt is relatively small, since Cj+1,0 = Cj on the left-hand side of Eq. (7) is infinitesimal. As shown in Fig. 1, the simulation tends to have large head changes and cut through the peak of the bell-shaped C curves. For instance, when Dt = 8.2  10
5 days, I = 15 cm/d, the calculated head hj+1,0+1 is 10
2 cm. Apparently, the change of head is large and the convergence criterion is not fulfilled. When k = 1 (second iteration), the corresponding C j+1,1 is still very small (10
4 cm
1 in Fig. 1) while K j+1,1 increases significantly (by five orders of magnitude in Fig. 1). A large K and relatively large hydraulic gradient lead to a large qj+1,1+1 and the calculated head hj+1,1+1 at the surface node decreases dramatically (see the h value at k = 2). The effective saturation and the K value for the first node will shift between very small and large values between adjacent iteration levels. In contrast, the C value is kept very small at all iteration levels. As a result, the calculated head deviates from the solution at time step j + 1 (hj+1,⁄) more and more significantly during the iteration. Thus, the divergence of this infiltration problem is induced by the nonlinearities of both soil water capacity and hydraulic conductivity curves, especially the nonmonotonicity of C. The problem will deteriorate rapidly if the soil saturation status changes during the iteration, as indicated in Fig. 1 at k = 3. On the other hand, this problem will be alleviated if a less nonlinear soil hydraulic model (e.g., BrooksCorey model (Brooks and Corey, 1964)) is used to describe the soil retention curve and the unsaturated conductivity curve. The above analysis is based on a flux upper boundary condition. When a prescribed head boundary condition is used (e.g., h = 0 used by Huang et al. (1996)), the divergence problems are less sev-

2.4.2. Influencing factors for the numerical difficulties The size of time step can control the head change during one iteration. For specified flow condition, there is an allowed maximum time step size Dtm that can maintain the convergence of the current iteration. Obviously, Dtm will change with grid size, soil type, hydraulic functions, initial and boundary conditions. In the following text, we will analyze the sensitivities of Dtm on these influencing factors. While all the other factors remain unchanged as their background values, different values are assigned to the studied influencing factor, and corresponding Dtm changes are recorded. The background problem is infiltration (constant rate I = 15 cm/d) into homogeneous loam soil (Table 1), with initially dry condition (h = -10,000 cm throughout the profile). The hydraulic parameters of loam soil are selected as the background parameters for the sensitivity analysis since they are close to the medians of the van Genuchten parameters for all soils. Additional sensitivity analyses using other mean parameters not shown here indicate that other background parameters do not affect the conclusions here. The 100-cm column is uniformly discretized into 10  10 cm elements. Fig. 2 reveals the influences of initial pressure head (Fig. 2(a)), mesh size (Fig. 2(b)), infiltration rate (Fig. 2(c)), and hydraulic parameters (Fig. 2(d–h)) on the allowed maximum time step size Dtm that can maintain the convergence of CI method. It is found that Dtm increases linearly with the size of the mesh (Fig. 2(a)). In other words, the problem with fine mesh is more difficult to ensure convergence, which is unfavorable since we need a fine mesh to reduce the spatial discretization error, especially at the soil surface where the hydraulic gradient is high (Vogel and Ippisch, 2008). The influence of initial condition (h0) on the required time step size has been investigated by Hills et al. (1989) and Zha et al. (2013b), and the results here (Fig. 2(b)) are consistent with their findings. Drier soil is more difficult to converge and requires smaller Dtm. It is also intuitive to find that Dtm decreases with the increase of infiltration rate I (Fig. 2(c)). Moreover, we also investigate the impacts of the hydraulic parameters. As shown in Fig. 2(d, e, f), Dtm decreases with the increases of parameters hr, a, n. Since hr is inversely proportional to the value of moisture capacity C, a large hr means that there is smaller space for the storage of infiltrated water and the surface node will likely become saturated for the first iteration. a is a nonlinear parameter that affects the moisture capacity C. Large a value leads to small C when the soil is dry. Thus, the infiltration problem tends to diverge for large a value. n is a nonlinear parameter that affects both the moisture capacity C and the unsaturated conductivity K. Large n

60

Y. Zha et al. / Journal of Hydrology 551 (2017) 56–69

Fig. 2. The influences of initial pressure head h0, infiltration rate I, mesh size Dz, and hydraulic parameters on the allowed maximum time step size Dtm for the convergence of Celia Picard iteration method. The background problem is the constant infiltration rate with I = 15 cm/d into soil with length of 600 cm. Bottom boundary is no-flux condition. The soil profile is discretized uniformly with Dz = 10 cm. Initial head is h0 =
10,000 cm over the profile. The background parameters are a = 0.036 cm
1, n = 1.56, hr = 0.078, hr = 0.43, and Ks = 24.96 cm/d. The influence of investigated factor is then obtained by changing the single factor and obtaining corresponding Dtm.

leads to very steep curve for the unsaturated conductivity, which means that the K value contrast will be enormous for the wet and dry situations, and this will induce the divergence of the simulation. Dtm increases with the increase of parameter hs (Fig. 2(g)), since hs is linearly proportional to C. The result also indicates that the saturated conductivity Ks does not have significant impact on the allowed Dtm (Fig. 2(h)) if the Ks change is not enormous. This is because the soil is relatively dry and has very low conductivity K initially. A linear factor (Ks) does not alter the steep shape of Kh curve, and the iterations will still jump between wet and dry statuses. In summary, it is more difficult to simulate the infiltration into dry, coarse-textured soil with a fine discretization and a high

infiltration flux intensity. Unfortunately, fine discretization is often required for accuracy purpose, and infiltration into the dry soil is very common due to atmospheric boundary condition.

3. Modified iteration algorithm As shown in the previous section, the simulation divergence is due to the overestimation of the head at the first iteration (hj+1,1  hj+1,⁄). The reason is that we estimate the change of head based on a very small moisture capacity (Cj+1,k) used in Picard iteration, while the average moisture capacity during time step j

61

Y. Zha et al. / Journal of Hydrology 551 (2017) 56–69

should be

1 hjþ1;kþ1
h j

R hjþ1;kþ1 hj

CðhÞdh = (hj+1,k+1
hj)/(hj+1,k+1
hj), which

is much larger than Cj+1,k. To prevent divergence, we need to (1) decide whether the iteration tends to fail, and if so, (2) make a correction to the guessed head values. Suppose at a certain time step j and iteration step k, we have obtained the intermediate head solu~ jþ1;kþ1 (i = 1, 2, . . ., N, N is the total number of nodes) by solvtion h i

ing the linear equation set. A modification to the CI scheme is proposed in Algorithm 1. For the first step, two user-specified parameters, hdry and hwet, are introduced to judge if the calculated head using CI method ~jþ1;kþ1 ) needs to be reevaluated or not (Fig. 3). There are three sce(h

jþ1;k

narios in this figure: (1) The soil is not very dry (hi > hdry ) and thus the initial C is not small enough to produce enormous head change during iterations; (2) The soil is initially very dry jþ1;k

6 hdry ) but the time step size is adequately small to get a (hi proper hj+1,k+1; (3) The soil is initially dry and the calculated head ~ jþ1;kþ1 P hwet ) is very large. Among them, scenario (3) needs a (h i

~ jþ1;kþ1 based on Cj+1,k is much recalculation since the guessed head h i larger than the convergent solution hj+1,⁄ at time step j + 1. In this study, hdry is chosen as the head when the soil is dry (on the left side of the hpeak in Fig. 3, where hpeak is the head when C is maximum) and the corresponding C is equal to 10
3 m
1 (or 10
5 cm
1). Similarly, hwet is chosen as the wet head (on the right side of hpeak in Fig. 3) when the corresponding moisture capacity value is equal to half of the maximum one. The hdry and hwet values related to different soils used in the numerical cases are listed in Table 1. The calculated hdry ranges from
3500 cm to
210 cm, with larger values for sandy soils. In contrast, hwet values are all close to zero (
2.6 cm to 0 cm). For other soil water retention curves (e.g., Brooks-Corey model) which do not have bell-shaped C, hwet can be simply set as zero. In summary, in our proposed algorithm the head values before jþ1;k ~jþ1;kþ1 ) the iteration are examined. If the original ) and after (h (h i

i

jþ1;k

is smaller than a user-specified threshold hdry, and the ~ jþ1;kþ1 is larger than a user-specified threshold h , updated head h wet i i.e.,

head hi

jþ1;k

hi

~jþ1;kþ1 P hwet ðM Þ; 6 hdry ðM i Þ and h i i

ð8Þ

the risk of numerical divergence for this node is deemed to be very high. Here Mi is the soil type index for node i. Under this situation, ~jþ1;kþ1 is firstly used to estimate the obtained intermediate head h i

: the soil moisture hjþ1;kþ1 i Fig. 3. The illustration of Eq. (8) and parameters hwet and hdry in MI algorithm. Scenario (1) has initially wet status and scenario (2) has an adequately small time step size. Scenario (3) satisfies Eq. (8) and recalculation is needed. Parameters hwet and hdry are heads when the moisture capacity values are 10
5 cm
1 and half of the maximum C value, respectively.

~jþ1;kþ1
hjþ1;k Þ hjþ1;kþ1 ¼ hjþ1;k þ C jþ1;k ðh i i i i i jþ1;kþ1

and the head hi

ð9Þ

at iteration level k + 1 is given via,

Table 2 Computational costs (cumulated iteration number) and convergence performances of MI and CI algorithms in the three numerical cases. Case

Case 1

Case 2

Case 3

Scenario

5000-1-10E-41 50000-1-E-4 5000-0.1-E-4 50000-0.1-E-4 50000-1-E-2 50000-0.1-E-2 Sand Loamy sand Sandy loam Loam Silt Silt loam Sandy clay loam Clay loam Silty clay loam Sandy clay Silty clay Clay /

Cumulated Iterations

Robustness

MI

CI

MI

CI

1748 1895 7174 22942 574 3443 17849 17415 20739 17347 17824 17242 20973 16170 16222 17989 15735 17675 22830

1747 1835 6768 21503 619 NA2 14343 18936 27432 18569 19075 18476 21733 16801 16794 20359 15852 17686 NA

Good4

Moderate3 Moderate Moderate Moderate Moderate Poor5 Poor Moderate Poor Moderate Moderate Moderate Moderate Good Good Moderate Good Good Poor

5000-1-10E-4 denotes initial head h0 =
5000 cm, grid size Dz = 1 cm and converging parameter eh = 10
4. Not available due to early termination of simulation using CI algorithm. 3 The simulation has no problem if proper initial time step size is carefully chosen. 4 The simulation has no problem without changing default time stepping parameters. 5 NaN values or saturated moisture profiles are obtained due to numerical instability; or the simulation cannot continue with several attempts of different time stepping parameters. 1 2

62

Y. Zha et al. / Journal of Hydrology 551 (2017) 56–69

jþ1;kþ1

hi

¼ hðhjþ1;kþ1 Þ i

ð10Þ

Algorithm 1. The proposed modified Picard iteration scheme (MI). At a particular time level j and iteration level k, solve the ~jþ1;kþ1 (i = 1, 2, . . ., N) from the linear equaupdated head values h i

tion set Ah = b based on Eq. (4). Do i=1 to N ~jþ1;kþ1 P hwet ðM Þ) Then 6 hdry ðM i Þ and h i i ~ jþ1;kþ1
hjþ1;k Þ hjþ1;kþ1 ¼ hjþ1;k þ C jþ1;k ðh

jþ1;k

If (hi

i jþ1;kþ1

hi

i

i

i

i

¼ hðhjþ1;kþ1 Þ i

Else jþ1;kþ1

hi Endif Enddo

~jþ1;kþ1 ¼h i

tion level k + 1. If Eq. (8) is not satisfied, CI updating scheme is ~ jþ1;kþ1 ? hj+1,k+1). As a result, MI algorithm only changes retained (h the guessed head at the first iteration to improve the convergence. Thus, it is expected to have similar accuracy, computational cost, and mass conservation performances as those of CI algorithm. Considering that the difference between MI and CI algorithms is not related to spatial discretization, it is reasonable to test the proposed method using only finite element model in the following numerical examples. The above procedure goes through all N nodes. Since usually only a minority of nodes (e.g., nodes at the soil surface) satisfy Eq. (8), most of the nodes still update the calculated head and soil moisture values according to Celia et al. (1990). From an implementation point of view, this MI algorithm only requires adding several lines to the original code using CI algorithm. 4. Numerical experiments

j+1,k+1

j+1,k+1

Eq. (9) indicates that we recalculate the new head h by h with the help of soil water retention curve and the intermediate ~ jþ1;kþ1 ) from solving the equation set is discarded. The varihead (h ables hj+1,k+1 and hj+1,k+1 are used to form the equation set for itera-

Based on different forms of Richards’ equations, Zha (2014) developed a series of finite element codes to solve 1-D, 2-D, or 3D variably saturated flow problems, and the codes were verified with analytical solutions and numerical solutions from other

Fig. 4. The evolution of time step size Dt versus time step level j for different scenarios (with different mesh sizes or initial conditions) simulated by CI and MI algorithms in case 1, which simulates infiltration into the layered soil with a constant rate of 2 cm/d, and initial head h0 =
5000 or
50,000 cm.

Y. Zha et al. / Journal of Hydrology 551 (2017) 56–69

numerical models, such as HYDRUS. The CI algorithm was implemented in the codes to ensure mass conservation. Starting from CI algorithm, we implement the proposed MI algorithm, as described by Eqs. (8–10). The usefulness of the proposed approach is demonstrated by comparing the robustness, efficiency, massconservation, and accuracy of the MI and CI algorithms. 4.1. 1-D infiltration problem The profile consists of five layers that alternate between two soil types, sand and loam (see Table 1), starting with the sand at the top. The total profile is 100 cm with each soil layer having an equal thickness of 20 cm. The simulation lasts for 5 days. The upper boundary is imposed with constant infiltration (I = 2 cm/d) and lower boundary is impermeable. We use two different initial heads (h0 =
5000 or
50,000 cm) and two grid sizes (Dz = 1 cm or 0.1 cm) to investigate their influences on the simulation performance. Moreover, for the scenarios with h0 =
50,000 cm, two different convergence tolerances (eh = 10
4 or 10
2) are used to investigate their impact on the mass conservation. Thus, there are a total of 6 scenarios (Table 2). This example is similar to Hills et al. (1989), except that the soil here is more challenging due to larger a value. The reference solutions are obtained by hform Richards’ equation on a very fine grid (Dz = 0.01 cm), since this type of model is more efficient as long as soil profile is unsaturated and the discontinuities between soil layers are carefully handled (Zha et al., 2013b).

63

It is obvious that the finer discretization leads to the larger size of equation sets, and thus requires more CPU time per iteration. Moreover, as illustrated in the previous section, fine grid requires smaller initial time step for convergence. As shown in Fig. 3, the initial time step size for algorithm CI should be smaller with the drier initial condition to avoid divergence (10
6 d for h0 = -5000 cm and 10
8 for h0 = -50,000 cm when Dz = 1 cm in Fig. 3). As a result, CI requires approximately 100 times the CPU as soil become drier if uniform time stepping is adopted. Because the soil becomes wetter in the process of infiltration, the numerical model can increase the time step size as time progresses. Due to adaptive time stepping, the algorithms CI and MI almost have the same iteration numbers when h0 =
5000 cm and Dz = 1 cm (Table 2). Similarly, CI also requires smaller Dtm when Dz is smaller (as indicated by the initial Dt in Fig. 4). In contrast, the initial time step size of the modified algorithm (MI) is not sensitive to the initial condition or Dz. With dryer initial condition or finer grid size, CI requires more iterations than MI. For instance, when h0 = -50,000 cm and Dz = 1 cm, CI requires approximately 1.2 times iterations of MI (Table 2). Within the same scenario, it is found that CI and MI have indistinguishable head/moisture profiles (thus not shown here), which means that the accuracy is not sensitive to the selected algorithm. This is reasonable since MI only optimizes the first-cut guess of h value, and does not change the major computational procedures. On the other hand, the accuracy is related to the spatial discretization. Comparing to the reference solution obtained with very fine

Fig. 5. The simulated moisture and head profiles at t = 2 d or 5 d using MI algorithm in case 1, which simulates infiltration into the layered soil with a constant rate of 2 cm/d, initial head h0 =
50,000 cm, and mesh size Dz = 1 cm or 0.1 cm. The reference solutions are obtained by h-form Richards’ equation on very fine grid (Dz = 0.01 cm).

64

Y. Zha et al. / Journal of Hydrology 551 (2017) 56–69

grid, the solution using coarse grid tends to have deeper infiltration wetting fronts (Fig. 5). The mass-conservation performances of MI and CI algorithms are also investigated in this case. Fig. 6 presents the mass balance errors of MI and CI algorithms versus time with different mesh sizes (Dz = 1 cm or 0.1 cm) or converging parameter (eh = 10
4 or 10
2). The mass balance error is given by the storage change of the soil moisture profile subtracting the amount of infiltration, i.e.,

EðtÞ ¼ WðtÞ
Wð0Þ
It

Even though MI and CI gain indistinguishable results and MI is only marginally superior in computational cost and mass conservation, the numerical stability of these two algorithms are quite different (Table 2). CI fails in the first attempt of simulation when the initial time step is set to be larger than the required time step size (Dtm) for convergence. For instance, with relatively dry initial condition (h0 = -50,000 cm) and fine discretization (Dz = 0.1 cm), it is found that the initial time step size has to be set less than 10
10 days for convergence, which is much smaller than the default value of time step size for common variably saturated flow models. For new modeler without considerable modeling experience, it is quite difficult to avoid divergence problem in this case. In contrast, the required time step for MI is not sensitive to initial condition or spatial discretization, thus it has more stable performance.

ð11Þ

RL where W(t) = 0 hðz; tÞdz [L] is the total water at time t over the soil profile with z2[0, L]. As shown in Fig. 6(a), when Dz = 1 cm and eh = 10
4, the mass balance errors of MI and CI are negligible (<0.00015 cm) compared to the infiltration amount or storage change in the profile. Using a loose converging parameter (eh = 10
2) leads to a slight increase of mass error (Fig. 6(b)), which can be well explained by Eqs. (4) and (5). When the mesh is refined (Dz = 0.1 cm), the mass balance error of CI is larger than that of MI (Fig. 6(c)). This could be attributed to the fact that fine grid is more difficult to converge with CI algorithm. When eh is increased to 10
2, the mass balance error of MI does not show significant change, while CI exhibits enormous mass error. By manually checking, we find that CI obtains a fully-saturated profile after time step j = 4, which is physically incorrect and thus produce large mass balance error. Due to better converging property, MI algorithm is superior to CI regarding mass conservation performance.

4.2. 1-D problem with atmospheric boundary condition Case 2 is designed to examine the CI and MI algorithms when atmospheric boundary conditions with alternating dryingwetting processes are imposed to the soil surface. Different from case 1, the dynamic atmospheric condition in case 2 induces abrupt change of convergence property using CI algorithm. The dry soil is subject to upward flux (i.e., only evaporation) during dry spells, and the solution is easy to converge. Once the upper boundary switches to infiltration, the convergence property deteriorates rapidly. However, heuristic time stepping techniques use the iteration number of previous time step to adjust the current

(a)

(b)

−5

0

5

10

15

εθ=0.01 Δz=1 cm

−10

−10

−5

0

5

10

15

εθ=0.0001 Δz=1 cm

1

2

3

4

5

(c)

1

2

3

4

5

2

3

4

5

(d) εθ=0.01 Δz=0.1 cm

150

εθ=0.0001 Δz=0.1 cm

150

0 200

200

0

0

0

50

50

100

100

Mass balance error (×10−5cm)

CI algorithm 20

20

MI algorithm

0

1

2

3

4

5

0

1

t (d) Fig. 6. The mass balance errors of CI and MI algorithms versus time in case 1 with different mesh sizes (Dz = 1 cm or 0.1 cm) or converging parameter (eh = 10
4 or 10
2). The mass balance error is given by the storage change of the soil moisture profile subtracting the amount of infiltration.

65

Y. Zha et al. / Journal of Hydrology 551 (2017) 56–69

The profile is 300 cm in length and consists of homogeneous soil. A total of 12 soils from Carsel and Parrish (1988) are adopted in this case (Table 1) to investigate the influence of soil texture on the simulation performance. The initial head profile is uniform and h0 =
5000 cm, and free drainage condition is specified for the bottom boundary. The potential evaporation and precipitation are displayed in Fig. 7, with a total simulation time of 365 days. Since the choices of grid size have been investigated in case 1, this case uses 101 nodes, with Dz = 5 cm at the bottom, and it is refined gradually to 0.1 cm on the top. The evolutions of time step size Dt in the simulations with 12 different soils using CI and MI algorithms are presented in Fig. 8. As expected, the time step size Dt using CI algorithm shows abrupt changes for the majority of the 12 soils, indicating that recalculation procedure (see Section 2.3) is triggered after simulation divergences, when water infiltrates into dry soil after a long period of evaporation. Among these scenarios, the simulation time step of sand using CI is initially very small and then it becomes normal (Fig. 8(a)). However, it is found that the simulated profile is wrong after the divergence (fully saturated profile, not shown here). On the other hand, the simulation of sandy loam soil is terminated abnormally when an infiltration process begins at t = 116 days (Fig. 8(c)). The simulations of (b) loamy sand, (d) loam, (e) silt, (f) silt loam, (g) sandy clay loam, and (j) sandy clay using CI method are successful and the calculated soil moistures are reasonable

t (d) Fig. 7. The potential evaporation and precipitation for case 2, where the homogenous soil column is imposed with atmospheric upper boundary and free drainage bottom boundary.

200

100

200

300

−5 −10 300

300

100

200

300

300

0

200

−5 0

100

200

300

(h) Clay loam 0

100

200

300

−5

−5 0

100

−10 (g) Sandy clay loam

−10

−15

(j) Sandy clay

0

0

200

−10 0

200

−5 100

0

0

−5

−5 −10 −15

(i) Silty clay loam

100

−10 (f) Silt loam

−15

300

0

0

200

−15

300

−10 100

0

0

0

0 −5 −10 100

−5

0 −5 −10

(e) Silt

−15

log10Δt (d)

0

(d) Loam

−10

300

0

200

(c) Sandy loam

−15

100

(b) Loamy sand

0

0

MI algorithm

(k) Silty clay

−15

−15

(a) Sand

−15

−10

−10

−5

−5

0

0

CI algorithm

−15

300

−15

200

0

100

200

300

(i) Clay

−15

10 100

15

0

Precipitation (cm)

5

1.5 1.0 0.5 0.0

Evaporation (cm)

2.0

0

time step size (see Section 2.3), which thus cannot accommodate the abrupt deterioration of convergence performance. In order to save the simulation from termination, CI has to recalculate this step with a much smaller Dt (smaller than the allowed maximum time step size Dtm in Section 2.4) once the model diverges.

0

100

200

300

t (d) Fig. 8. The evolutions of time step size Dt versus time t simulated by CI and MI algorithms in case 2, where the 12 homogeneous soil columns are imposed with atmospheric upper boundary and free drainage bottom boundary.

66

Y. Zha et al. / Journal of Hydrology 551 (2017) 56–69

and 0.2 days (Fig. 11). In contrast, the simulations using CI algorithm always terminate before the end of simulation time or obtain abnormal results (NaN or saturated moisture values) after several attempts have been tried (see Table 2). The divergence problem for simulating infiltration flow into the dry soil in case 3 is severer due to the complexity of the 2-D flow, as indicated by the soil moisture and corresponding head maps at t = 198 and 364 days in Fig. 12 and the head animation in supporting material. For instance, the wetting front at t = 198 days is deeper at x = 400– 500 cm, possibly due to the more-permeable soil at the top and less-permeable soil below the wetting front. As a result, the streamlines at the front indicate that the flow moves horizontally due to the barrier below the wetting front. Eventually, the wetting front arrives at the water table and the flow is mainly vertical in unsaturated zone towards groundwater and horizontal towards discharge boundaries in the saturated zone.

500

sandy loam

200 100

sandy clay loam

0

sand

0

200

400

600

800

x (cm) Fig. 10. The soil facies model for case 3 with heterogeneous soil classes. The 2-D domain consists of an 800 by 650 cm rectangular box. The top is subject to the atmospheric boundary, and the water table at the left and right sides is maintained at z = 200 cm.

MI algorithm 0.30

loam

300

Previous cases are limited to 1-D flow problems. Case 3 is designed to examine the CI and MI algorithms for simulating 2-D flow. The 2-D domain consists of a 800 by 650 cm rectangular box, filling with heterogeneous soil with four different facies, as shown in Fig. 10. The top soil is subject to the same atmospheric boundary condition used in case 2. Initially the soil has established hydrostatic equilibrium with the water table at 200 cm from the bottom. The water tables at the left and the right sides are maintained as 200 cm (i.e., constant head boundary). The total simulation time is 365 days. The 2-D domain is discretized into 90  49 meshes (91  50 nodes). The mesh in the horizontal direction is uniform with Dx = 8 cm. The mesh is gradually refined from bottom (Dz = 25 cm) to the top (Dz = 1 cm) since the hydraulic gradient induced by infiltration is higher at the surface. Again, the MI algorithm successfully simulates the soil water dynamics for one year. The time step size Dt is between 7  10
5

z (cm)

4.3. 2-D infiltration problem in heterogeneous soil

400

600

(by manual check), although they still suffer from intensive recalculations. The remaining simulations using CI are readily completed without demanding recalculations for fine-textured soils, including (h) clay loam, (i) silty clay loam, (k) silty clay, and (i) clay, where the values of parameter a are all smaller than 0.02 cm
1. In contrast, the evolution of Dt using MI algorithm is not sensitive to the soil textures, and the recalculation procedure is rarely provoked. The recalculation requires additional coding effort, and the iterations in divergence steps also cost more computational time. Sometimes the simulations are completed with incorrect moisture profiles (e.g., fully saturated or showing NaN (not a number) values) because CI is prone to produce extremely large head changes. Under this situation, the recalculation procedure using smaller time step size after divergence is not triggered and the final results require carefully manual checks. Thus, MI has much better convergence performance and it generally requires less iterations than CI (Table 2) even though the two algorithms achieve indistinguishable results (Fig. (Fig. 9)).

CI algorithm

0.30 0.00

0.15

(b) depth of 10 cm

(c) depth of 30 cm

0.15 0.00

θ

0.30 0.00

0.15

(a) depth of 0 cm

0

100

200

300

365

Time (d) Fig. 9. The simulated soil moistures at depths of 0, 10, and 30 cm versus time by CI and MI algorithms in case 2 using loamy sand soil.

67

Δt (d)

1e−05

1e−04

1e−03

1e−02

1e−01

Y. Zha et al. / Journal of Hydrology 551 (2017) 56–69

0

2000

4000

6000 8000 Time step level j

10000

12000

Fig. 11. The evolutions of time step size Dt versus time step level j simulated by MI algorithm in case 3.

600

600

θ 0.42 0.38 0.34 0.3 0.26 0.22 0.18 0.14 0.1 0.06

400 300 200 100 0

500

Z (cm)

Z (cm)

500

200

400

600

400 300 200 100

(a) t= 0.198E+03 Days 0

θ 0.42 0.38 0.34 0.3 0.26 0.22 0.18 0.14 0.1 0.06

0

800

(b) t= 0.364E+03 Days 0

200

X (cm)

400

600

600

600 h (cm)

400 300 200 100

(c) t= 0.198E+03 Days 0

200

400

600

-10 -20 -30 -40 -50 -60 -70 -80 -90 -100 -110 -120 -130 -140 -150

800

X (cm)

h (cm)

500

Z (cm)

Z (cm)

500

0

800

X (cm)

400 300 200 100 0

(d) t= 0.364E+03 Days 0

200

400

600

-10 -20 -30 -40 -50 -60 -70 -80 -90 -100 -110 -120 -130 -140 -150

800

X (cm)

Fig. 12. The soil moisture and head maps at t = 198 d and 364 d simulated by MI algorithm in case 3. The black lines with arrows are streamlines.

5. Summary and conclusions It is well-known that the numerical solution of Richards’ equation is time-consuming and often confront divergence issue for simulating infiltration into the soils. This is quite frustrating for modelers who want to obtain accurate soil water dynamics. As a matter of fact, although variably saturated models have been developed for regional-scale problems (Twarakavi et al., 2008), their applications are rare due to model stability and other issues

such as lack of hydraulic parameter information (Niswonger and Prudic, 2009). Nevertheless, there is no doubt that as the advance of computer science, numerical techniques, as well as field soil moisture measurement methods, the Richards’ equation models are preferred to capture the accurate soil water dynamics than water balance models. In order to alleviate the model stability issue, we first analyze the model divergence problem by monitoring the changes of hydraulic parameters (unsaturated hydraulic conductivity and

68

Y. Zha et al. / Journal of Hydrology 551 (2017) 56–69

moisture capacity) with heads during the first three iterations. To quantitatively analyze the factors that influence the model stability, we investigate the required maximum time step to maintain simulation convergence for a problem simulating infiltration into the dry soil. After revealing the reasons for the divergence issue, a slight modification of the state variables (i.e., pressure head and soil moisture content) updating procedure in Celia Picard iteration method is proposed. We conduct a 1-D infiltration problem and more realistic 1-D and 2-D cases with evaporation and infiltration (atmospheric boundary) in the homogeneous or heterogeneous soils. The results demonstrate the effectiveness of the modified iteration algorithm. The following conclusions are drawn from our analysis and numerical experiments: (1) The nonlinearities in the soil moisture capacity and unsaturated hydraulic conductivity contribute to the divergence issue for infiltration into dry soil. The extremely small soil moisture capacities during very dry and wet statuses are sensitive to any small amount of inflow/outflow and lead to enormous head changes in the iterations. In the meantime, the unsaturated conductivity changes rapidly between dry and wet and thus leads to unstable outflow during the iterations. (2) The maximum allowed time step Dtm for convergence in Celia Picard iteration is negatively correlated to the values of van Genuchten parameters hr, a and n, but it is positively correlated to parameters hs, initial head, and spatial discretization size Dz. Thus, the divergence issue is severer in simulating infiltration into very dry, sandy soil with fine spatial discretization. (3) Compared to the Celia Picard iteration method, the proposed modification can significantly enhance the model robustness, and it can marginally reduce the computational time and mass balance error, without losing any accuracy. The superiority of the proposed algorithm is more significant in coarse-textured soil.

Acknowledgments The authors acknowledge the supports from National Natural Science Foundation of China (No. 51609173, 51479144, 51522904), the open-end fund provided by Key Laboratory for Groundwater and Ecology in Arid and Semi-Arid Areas, CGS (No. KLGEAS201601), and the Special Fund for Public Industry Research from Ministry of Land and Resources of China (201511047). The valuable suggestions from two anonymous reviewers are appreciated. We thank Deqiang Mao at New Mexico Tech for proofreading the manuscript. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.jhydrol.2017.05. 053. References Brooks, R.H., Corey, T., 1964. Hydraulic Properties of Porous Media. Colorado State University, Fort Collins. Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water retention characteristics. Water Resour. Res. 24, 755–769. Caviedes-Voullième, D., Garcı’a-Navarro, P., Murillo, J., 2013. Verification, conservation, stability and efficiency of a finite volume method for the 1D Richards equation. J. Hydrol. 480, 69–84. http://dx.doi.org/10.1016/j. jhydrol.2012.12.008. Celia, M., Bouloutas, E., Zarba, R., 1990. A general mass-conservative numerical solution for the unsaturated flow equation. Water Resour. Res. 26, 1483–1496.

Crevoisier, D., Chanzy, A., Voltz, M., 2009. Evaluation of the Ross fast solution of Richards’ equation in unfavourable conditions for standard finite element methods. Adv. Water Resour. 32, 936–947. http://dx.doi.org/10.1016/j. advwatres.2009.03.008. D’Haese, C.M.F., Putti, M., Paniconi, C., Verhoest, N.E.C., 2007. Assessment of adaptive and heuristic time stepping for variably saturated flow. Int. J. Numer. Meth. Fluids 53, 1173–1193. http://dx.doi.org/10.1002/fld.1369. Diersch, H.J.G., Perrochet, P., 1999. On the primary variable switching technique for simulating unsaturated-saturated flows. Adv. Water Resour. 23, 271–301. http://dx.doi.org/10.1016/S0309-1708(98)00057-8. Dogan, A., Motz, L.H., Asce, M., 2005. Saturated-unsaturated 3D groundwater model. I: Development. J. Hydrol. Eng. 10, 492–504. Forsyth, P., Wu, Y., Pruess, K., 1995. Robust numerical methods for saturatedunsaturated flow with dry initial conditions in heterogeneous media. Adv. Water Resour. 18, 25–38. Hao, X., Zhang, R., Kravchenko, A., 2005. A mass-conservative switching method for simulating saturated–unsaturated flow. J. Hydrol. 311, 254–265. http://dx.doi. org/10.1016/j.jhydrol.2005.01.019. Hassane Maina, F., Ackerer, P., 2016. Ross scheme, Newton-Raphson iterative methods and time-stepping strategies for solving the mixed-form of Richards’ equation. Hydrol. Earth Syst. Sci. Discuss. 1–48. http://dx.doi.org/10.5194/hess2016-622. Hills, R., Porro, I., Hudson, D.B., Wierenga, P.J., 1989. Modeling one-dimensional infiltration into very dry soils: 1. Model development and evaluation. Water Resour. Res. 25, 1259–1269. Huang, K., Mohanty, B., Genuchten, M..Van., 1996. A new convergence criterion for the modified Picard iteration method to solve the variably saturated flow equation. J. Hydrol. 178, 69–91. Kavetski, D., 2002. Noniterative time stepping schemes with adaptive truncation error control for the solution of Richards equation. Water Resour. Res. 38. http://dx.doi.org/10.1029/2001WR000720 (29-(1–10)). Kavetski, D., Binning, P., Sloan, S.W., 2001. Adaptive time stepping and error control in a mass conservative numerical solution of the mixed form of Richards equation. Adv. Water Resour. 24, 595–605. http://dx.doi.org/10.1016/S03091708(00)00076-2. Kirkland, M.R., Hills, R.G., Wierenga, P.J., 1992. Algorithms for solving Richards equation for variably saturated soils. Water Resour. Res. 28, 2049–2058. Kosugi, K., 2008. Comparison of three methods for discretizing the storage term of the Richards equation. Vadose Zo. J. 7, 957–965. http://dx.doi.org/10.2136/ vzj2007.0178. Krabbenhøft, K., 2007. An alternative to primary variable switching in saturated– unsaturated flow computations. Adv. Water Resour. 30, 483–492. http://dx.doi. org/10.1016/j.advwatres.2006.04.009. Lai, W., Ogden, F.L., 2015. A mass-conservative finite volume predictor – corrector solution of the 1D Richards ’ equation. J. Hydrol. 523, 119–127. http://dx.doi. org/10.1016/j.jhydrol.2015.01.053. Mao, D., Wan, L., Yeh, T.-C.J., Lee, C.-H., Hsu, K.-C., Wen, J.-C., Lu, W., 2011. A revisit of drawdown behavior during pumping in unconfined aquifers. Water Resour. Res. 47, W05502. http://dx.doi.org/10.1029/2010WR009326. Matthews, C.J., Braddock, R.D., Sander, G.C., 2004. Modeling flow through a onedimensional multi-layered soil profile using the Method of Lines. Environ. Model. Assess. 9, 103–113. Menziani, M., Pugnaghi, S., Vincenzi, S., 2007. Analytical solutions of the linearized Richards equation for discrete arbitrary initial and boundary conditions. J. Hydrol. 332, 214–225. http://dx.doi.org/10.1016/j.jhydrol.2006.06.030. Miller, C., Abhishek, C., Farthing, M., 2006. A spatially and temporally adaptive solution of Richards’ equation. Adv. Water Resour. 29, 525–545. http://dx.doi. org/10.1016/j.advwatres.2005.06.008. Niswonger, R.G., Prudic, D.E., 2009. Comment on ‘‘Evaluating Interactions between Groundwater and Vadose Zone Using the HYDRUS-Based Flow Package for MODFLOW” by Navin Kumar C. Twarakavi, Jirka Simunek, and Sophia Seo. Vadose Zo. J. 8, 818–819. http://dx.doi.org/10.2136/vzj2008.0155. Rathfelder, K., Abriola, L.M., Arbor, A., 1994. Mass conservative numerical solutions of the head based Richards equation. Water Resour. Res. 30, 2579–2586. Ross, P.J., 2003. Modeling soil water and solute transport—fast, simplified numerical solutions. Agron. J. 95, 1352–1361. Sadegh Zadeh, K., 2011. A mass-conservative switching algorithm for modeling fluid flow in variably saturated porous media. J. Comput. Phys. 230, 664–679. http:// dx.doi.org/10.1016/j.jcp.2010.10.011. Simunek, J., Sejna, M., Saito, H., Sakai, M., van Genuchten, M.T., 2009. The HYDRUS1D Software Package for Simulating the One-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media. Department of Environmental Sciences, University of California Riverside, Riverside, CA. Srivastava, R., Yeh, T.-C.J., 1991. Analytical solutions for one-dimensional, transient infiltration toward the water table in homogeneous and layered soils. Water Resour. Res. 27, 753–762. Srivastava, R., Yeh, T.-C.J., 1992. A three-dimensional numerical model for water flow and transport of chemically reactive solute through porous media under variably saturated conditions. Adv. Water Resour. 15, 275–287. http://dx.doi. org/10.1016/0309-1708(92)90014-S. Tracy, F.T., 2006. Clean two- and three-dimensional analytical solutions of Richards’ equation for testing numerical solvers. Water Resour. Res. 42, 1–11. http://dx. doi.org/10.1029/2005WR004638. Twarakavi, N.K.C., Šimu˚nek, J., Seo, S., 2008. Evaluating interactions between groundwater and vadose zone using the HYDRUS-based flow package for MODFLOW. Vadose Zo. J. 7, 757–768. http://dx.doi.org/10.2136/vzj2007.0082.

Y. Zha et al. / Journal of Hydrology 551 (2017) 56–69 van Dam, J.C., Feddes, R.A., 2000. Numerical simulation of infiltration, evaporation and shallow groundwater levels with the Richards equation. J. Hydrol. 233, 72– 85. van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898. Varado, N., Braud, I., Ross, P.J., 2006a. Development and assessment of an efficient vadose zone module solving the 1D Richards’ equation and including root extraction by plants. J. Hydrol. 323, 258–275. http://dx.doi.org/10.1016/j. jhydrol.2005.09.015. Varado, N., Braud, I., Ross, P.J., Haverkamp, R., 2006b. Assessment of an efficient numerical solution of the 1D Richards’ equation on bare soil. J. Hydrol. 323, 244–257. http://dx.doi.org/10.1016/j.jhydrol.2005.07.052. Vogel, H.-J., Ippisch, O., 2008. Estimation of a critical spatial discretization limit for solving richards’ equation at large scales. Vadose Zo. J. 7, 112–114. http://dx.doi. org/10.2136/vzj2006.0182.

69

Zha, Y., 2014. Research on Cost-Effective Algorithm for Unsaturated-Saturated Flow and Its Application. Wuhan University. Zha, Y., Shi, L., Ye, M., Yang, J., 2013. A generalized Ross method for two- and threedimensional variably saturated flow. Adv. Water Resour. 54, 67–77. http://dx. doi.org/10.1016/j.advwatres.2013.01.002. Zha, Y., Yang, J., Shi, L., Song, X., 2013. Simulating one-dimensional unsaturated flow in heterogeneous soils with water content-based Richards equation. Vadose Zo. J. 12, 1–12. http://dx.doi.org/10.2136/vzj2012.0109. Zha, Y., Tso, M.C.-H., Shi, L., Yang, J., 2016. Comparison of noniterative algorithms based on different forms of Richards’ equation. Environ. Model. Assess. 21, 357– 370. http://dx.doi.org/10.1007/s10666-015-9467-1. Zhang, Z., Wang, W., Chen, L., Zhao, Y., An, K., Zhang, L., Liu, H., 2015. Finite analytic method for solving the unsaturated flow equation. Vadose Zo. J. 14, 1–10. http:// dx.doi.org/10.2136/vzj2014.06.0073.