Abstract

Key Words

Layered soils, discontinuity, water content, numerical solution, unsaturated flow.

Introduction

The study of water flow through layered soils is an important environmental problem, which has implications in various fields such as contaminant transport, agricultural science and waste disposal (Heilig et al. 2003; Colominas et al. 2002). Mathematical modelling has been used extensively as a tool to understand water flow processes through multiple soil layers (Srivastava and Yeh 1991; Warrick et al. 1997). Such models are generally based on the hydraulic pressure (h) form of Richards’ Equation since h is continuous across the interface between two soils. The main disadvantage of the h form of Richards Equation is that special consideration must be given to ensure that mass is conserved through the full wetting range especially for infiltration into dry soils (e.g. chord sloping approximation (Ratherfelder and Abriola 1994)). In contrast, the water content (θ) form of Richards’ Equation is mass conserving and accurate for infiltration into dry soils. However, the main disadvantages are that θ cannot be used to model saturated flow, numerical difficulties can arise for near-saturated conditions and θ is discontinuous across soil layer boundaries (Diersch and Perrochet 1999; Celia et al. 1990). This paper will address the latter problem.

Hills et al. (1989) were the first to attempt a θ-based model for infiltration into a layered soil profile. The solution was based on a Crank-Nicolson finite differencing scheme, which used a spatial grid that situated the soil layer interface halfway between two nodes. Hills et al. (1989) then used a transformation of Richards’ Equation, based on a step function, to specify a jump condition between two nodes near the interface. As commented by Hills et al. (1989), this approach will cause a small error since hydraulic properties are evaluated a half step length away from the interface. This error was further analysed by Romano et al. (1998) who concluded that the Hills et al. (1989) approach could produce an error at the interface as high as 4% when compared to the analytical solution of Srivastava and Yeh (1991). Note that Romano et al. (1998) developed an iterative version of the Hills et al. (1989) solution showing excellent agreement with the analytical solution but did not extend the solution for the q-form of Richards’ Equation. Additionally, finite element methods have been used to model water flow through heterogeneous soils; however, these models rely on averaging techniques to ensure continuity of flux across element edges, which can result in water balance errors (Diersch and Perrochet 1999; Diaw et al. 2001).

Recently, Matthews et al. (2004) developed a solution for θ-form of Richards’ Equation to simulate one-dimensional water flow through a multi-layered soil using the Method of Lines (MoL). The MoL, developed in Matlab, used a finite differencing scheme from Schiesser (1991) to descritise the spatial component of Richards’ Equation while keeping the time component continuous. This resulted in a system of Ordinary Differential Equations (ODEs) with each ODE representing the temporal gradient at a particular node within the spatial domain. For each test case in this study, the ODE system is integrated using a stiff ODE solver (ODE15s) from Matlab. The ODE15s integrator is a variable time stepping scheme based on Gears Method. As depicted in Figure 1, the spatial grid was constructed with a node situated on the interface boundary at z = z_c. A temporary node was also introduced on the interface allowing both soils to be represented at the boundary. This facilitates an approximation to the interface boundary conditions, which was used to construct an ODE at z = z_c giving an expression for at the interface. The MoL solution compared well against the analytical solution of Srivastava and Yeh (1991) and a numerical test case from Hills et al. (1989), but still relied on an extra approximation.

Figure 1. Grid system for 1D layered soil (Matthews et al. 2004)

This paper will use the MoL model of Matthews et al. (2004) and incorporate an iterative solution based on Newton’s Method that explicitly solves for θ at each soil layer interface. The solution will not rely on any additional approximations or assumptions at the interface except those already associated with the finite differencing scheme and the ODE integrator. It will be shown that the Newton’s function (f) is highly effective converging within 2 to 3 iterations. To ensure the convergence behaviour holds, the numerical solution will be tested on a fine over coarse soil and a coarse over fine soil with highly contrasting soil properties. Both soil profiles will be subjected to different surface flux intensities and the contrast between each soil type will be artificially controlled to extend and decrease the extent of the θ discontinuity. The new MoL model will also be compared against the model from Matthews et al. (2004), using an example from Hills et al. (1989), to explore advantages or disadvantages of using the iterative solution in the MoL. This paper will also ensure that the correct steady state is obtained by comparing the new MoL solution against the analytical solution of Warrick et al. (1997).

Governing equations

To simulate one-dimensional vertical flow through a layered soil profile as depicted in Figure 1, the problem can be stated mathematically in terms of θ as subject to the initial and boundary conditions

(1)

,

(2)

,

(3)

,

(4)

,

(5)

where k (= 1 or 2) is the soil layer counter, D (= K/C) is the soil-water diffusivity (L²T^-1), K is the hydraulic conductivity (LT^-1), is the water capacity (L^-1), q_s is a constant surface flux (LT^-1), z is the vertical spatial coordinate defined as positive downwards (L), t is time (T) and q_k is flux (LT^-1) defined as

.

(6)

In addition, the water retention curve, the hydraulic conductivity and hydraulic diffusivity for each soil will be represented by the van Genuchten (1981) equations as

,

(7)

,

(8)

,

(9)

respectively, where Θ is scaled water content, θ_r is the residual water content (L³L^-3), θ_s is the saturated water content (L³L^-3), K_s is the saturated hydraulic conductivity and n, m and α (L^-1) are curve fitting parameters related to soil texture. Note that equation (3) represents a constant flux boundary condition at the soil surface (z = 0) and the bottom boundary condition (equation 5) is a free drainage boundary condition where
q₂ = K₂(θ₂) at z = L.

Interface Boundary condition in terms of q

The interface boundary condition at z = z_c is based on the underlying physical principles of Richards’ Equation, that is, continuity of pressure and mass. This is represented in (4) by equating h and q at z = z_c. Both equations in (4) can be written explicitly in terms of θ by 1) rearranging equation (7) in both soils for h and substituting the resultant equation into the pressure equation and 2) substituting (6) into the flux equation, which gives respectively

,

(10)

,

(11)

where subscript 1 and 2 denote Soil 1 and 2 (Figure 1). Note that (10) represents the extent of the discontinuity in water content at the layer interface (z = z_c) since it is derived by equating h in both soils at
the interface. Since we have expressed (4) in terms of θ, it is now possible to derive an iterative scheme to solve for θ at the interface between the two soils. However, as a consequence of the discontinuity in θ, there are two values of θ at z = z_c i.e. θ₁ and θ₂. To implement the iterative scheme, a single solution variable must be chosen and we will arbitrarily choose θ₁. Note that θ₂ can be calculated from θ₁ through the relationships specified in (7) and (10).

Iterative Scheme using Newton’s Method

For this study, the iterative solution will be based on Newton’s Method, which involves finding the root of a function, that is, f = 0 using (12)

,

(12)

where b is the iteration counter, θ₁ is the solution variable, f is the Newton’s function and f ¢ is the derivative of f with respect to the solution variable θ₁. To generate a solution for θ₁ at the interface, we will solve (11) iteratively subject to the relationship in (10). Therefore, f can be written as

.

(13)

Within the MoL scheme, and at z = z_c are calculated using an upwinding and downwinding differencing equations (Matthews et al. 2004; Schiesser 1991) and will, therefore, be represented using Δ_u and Δ_d, respectively. Using this notation and (7) and (10), (13) can be expressed in discretised form as

.

(14)

Taking the derivative of (14) with respect to the solution variable θ₁ yields

.

(15)

Note that the functional notation in (15) is dropped and the derivatives , , are derived from (8), (9)and (10) in conjunction with (7). In summary, the complete iterative scheme, based on Newton’s method, is given by (7) to (15).

Implementation in the MoL

Fundamentally, the iterative solution converts the interface boundary conditions from a Neumann to a Dirichlet type boundary condition since a value for θ₁, and consequently θ₂ from (10), at z = z_c is obtained from (12) at each time step. In the MoL, the new θ₁ and θ₂ values are then imposed at node c in the water content vectors, θ₁ and θ₂, which contain θ at each node in Soil 1 and 2, respectively. Note that both θ₁ and θ₂ contain node c to represent both soils at the interface. The θ vectors are then used in subsequent calculations to provide a discretised form of (1) resulting in an ODE for each node in the system. For a Dirichlet boundary condition, it is recommended that at node c are set to zero to restrain the integrator from evaluating θ at the interface, thereby ensuring that the integrator does not perform any unnecessary work. For further details on MoL refer to Matthews et al. (2004). Note that the initial starting value of θ₁ for the iterative solution is taken from the previous time step and the convergence criteria for the iterative solution is set at .

Test Case 1 and 2: Explore the nature of the iterative solution

The first two test cases considered will comprise a fine over coarse soil (Test Case 1) and a coarse over fine soil (Test Case 2) taken from Hills et al. (1989). The soil parameters for these soils are provided in Table 1 and are highly contrasting. Both test cases will be used to explore the shape of the Newton’s function (f) and determine whether the convergence of the iterative scheme is affected by different soil stratifications and wetting events. In particular, the rate of convergence of the iterative scheme will be compared for different surface fluxes (q_s) while keeping the same initial water content distribution. As an additional test, the contrast between the two soils will be artificially controlled to explore the effects of decreasing and increasing the contrast between the two soils on the iterative solution. This will be achieved by defining a parameter ε so that the properties of the coarse soil will approach the properties of the fine soil as ε → 1. This relationship is given by and is applied to each soil parameter p in Table 1. Note that will be used to explore the effect of decreasing the contrast between the two soils where ε = 1 results in a homogenous soil profile. Similarly, ε < 0 will be used to explore the effect of increasing the soil contrast by increasing the coarseness of the coarse soil. Note that for both test cases, the depth of the soil profile is set at 40 cm with the interface at 20 cm, the initial condition corresponds to
h = -1000 cm and the spatial step size Δz = 0.25cm.

Table 1. Soil parameters for soils taken from Hills et al. (1989).

Soil type	Ks (cm/day)	θr	θs	α (cm-1)	n
Berino fine sand	541.0	0.0286	0.3658	0.0280	2.2390
Glendale clay	13.1	0.1060	0.4686	0.0104	1.3954

Test Case 3: Comparison against approximate MoL solution

The third test case will involve a comparison between the approximate MoL solution from Matthews et al. (2004) and the iterative MoL solution using an example from Hills et al. (1989). The Hills et al. (1989) problem involves infiltration into a very dry soil consisting of five layers that alternates between the two soil types given in Table 1 starting with the coarse soil (Berino fine sand). The total length of the soil profile is 100 cm with each soil layer having equal depth of 20 cm. The initial and boundary conditions for the problem are given by (2) to (4) with q_s = 2 cm/day and initial water content, θ_i1 and θ_i2, corresponding to
h = -1000 cm. For the interface boundary conditions in (4), k = 1, 2, 3 and 4 correspond to each interface at
z = 20 cm, 40 cm, 60 cm and 80 cm, respectively. For this problem, the bottom boundary condition (5) is replaced with the condition that θ₂ = θ_i2 at z = 100 cm. To compare the two MoL models, relative differemce (RD) will be calculated for each node within the system as , where θ_a and θ_e are water contents for a particular node from approximate and iterative solution, respectively. By comparing the two MoL models, we will explore whether there is any improvement in the accuracy and efficiency of the numerical solution by using the iterative solution as opposed to the approximate solution for the interface boundary condition. As an additional test, the steady state analytical solution of Warrick et al. (1997) will be applied to the Hills et al. (1989) example and compared to the numerical steady state solution from the MoL for the iterative case only. For the MoL solution, the criterion for steady state is set at . The relative error for the steady state case is calculated as above with θ_a and θ_e corresponding to MoL and analytical solution, respectively. Note that, for Test Case 3, the spatial step size was set at Δz = 0.1cm and all simulations were run on a Dell Desktop computer with a Pentium 4 processor and 2 Gigabytes of RAM.

Figure 2. Plot of depth vs θ for (a) Test Case 1 and (b) Test Case 2.

Results

Test Case 1 and 2

Figure 2 shows θ distribution over depth for Test Case 1 and 2 at 0.5 day intervals for the first 3 days for
q_s = 2 cm/day. At t = 3 days, both profiles were approaching steady state. On comparing Figures 2(a) and (b), the different flow behaviour through the different soil types and stratifications are evident. In both (a) and (b), the shape of the wetting front is much sharper in the coarse soil as compared to the fine soil where dispersion is apparent. Water had reached the interface by t = 1 day in both cases (Figure 2). For case 2 (Figure 2(b)), water had moved into the underlying fine soil penetrating to a depth of approximately 5 cm. For the fine over coarse soil (Figure 2(a)), water did not penetrate the underlying coarse soil to a depth of 5 cm until t = 1.5 days. This flow behaviour is known as the capillary barrier effect where water is held at the interface by capillary forces (Warrick et al. 1997). The holding capacity of the fine soil has also caused the top layer in Figure 2(a) to approach steady state faster with a small change in the profile from t = 2 days to
t = 3 days.

Table 2. Iteration profile for Test Case 1 and 2.

	Iterations	No. Function Eval.	% of Total	Cumulative %
Test Case 1	0	1017	38.0	38.0
	1	880	32.9	70.9
	2	570	21.3	92.2
	3	205	7.7	99.9
	>3	4	0.1	100
Test Case 2	0	1524	54.6	54.6
	1	954	34.2	88.8
	2	66	2.4	91.2
	3	245	8.8	100
	>3	0	0	100

Within the MoL, the iterative solution is implemented every time the discretised version of (1) is evaluated. These function evaluations are used within Matlab to measure the amount of work done by an ODE integrator. Table 2 shows the number of function evaluations that resulted in 0, 1, 2, 3 and greater than 3 iterations of the iterative solution. From Table 2, it is evident that the iterative solution converges very quickly for both test cases with approximately 90% of all function evaluations converging within 2 iterations. Also, the total number of function evaluations is higher for Test Case 2 (2789) as compared to Test Case 1 (2676). This extra work is caused by infiltration into the coarse soil having a stronger convective flow component. Interestingly, for Test Case 1, the majority of iterations greater than or equal to 3 were concentrated around the time period when the wetting front hit the interface. For Test Case 2, this was not the case where iterations greater than or equal to 3 were scattered throughout the entire simulation i.e. no apparent clustering.

Figure 3. Plot of f vs θ over time for (a) Test Case 1 and (b) Test Case 2. The circles highlight the root at f = 0.

The quick convergence of the iterative scheme is related to the shape of the Newtons function (f) over time but more importantly, the appropriateness of the initial guess to start the iterative scheme. Figure 3 shows the shape of the f for t = 1, 2 and 3 days for Test Case 1 and 2 with the solution at f = 0 highlighted by a circle. For both test cases, the shape of f is very similar, with a local maximum becoming more pronounced as time evolves, which is more predominant in Test Case 2. For most time intervals in Figure 3, the solution to the iterative scheme is situated on the downward arm of the function, which is near linear in nature. This behaviour is to a lesser extent for t = 1 day for both test cases in (a) and (b). Note that the initial guess for the iterative solution for all time intervals considered is situated very close to the solution and is not shown here since it is difficult to discern from the solution in Figure 3. Therefore, the key factors for the quick convergence is the close proximity of the initial guess to the actual solution and that f is a well behaved function near the solution throughout the entire simulation period for both cases. Note that after t = 3 days, there is very little change to the shape of f since the system is starting to approach steady state.

Figure 4. The effect of convergence on increasing surface flux (q_s) for (a) Test Case 1 and (b) Test Case 2.

Changes to the shape of f, or the positioning of the solution and the initial guess on f could be brought about by increasing the intensity of q_s or by increasing or decreasing the contrast between the two soil types. Figure 4 shows a plot of the total number of function evaluations (total work), number of function evaluations that had 0 to 2 iterations and greater than 2 iterations for both Test Case 1 and 2. Figure 4 (a) and (b) clearly shows that the convergence behaviour is preserved for increasing q_s with frequency of greater than 2 iterations staying reasonable constant. Figure 4(b) does show a slight increase in iterations greater than 2 but it is marginal compared to the total work done. Also, Figure 4 does show an increase in the total work done as q_s increases. This is to be expected since increasing q_s will result in sharper wetting fronts causing the integrator to take smaller time steps to keep within integrator tolerances.

Figure 5. The effect of convergence on increasing soil contrast (ε) for (a) Test Case 1 and (b) Test Case 2.

The effect of varying the contrast between the two soil properties via the ε parameter is shown in Figure 5 for Test Case 1 and 2. Figure 5 (a) and (b) clearly show that the number of function evaluations that result in iterations greater than 2 remain fairly constant for both test cases. This demonstrates that increasing or decreasing the contrast between the two soils does not have a large affect on the convergence behaviour of the iterative solution. Also, Figure 5 (a) and (b) show that the total number of function evaluations increases as ε decreases. This is to be expected since the model needs to handle a larger instantaneous jump in soil properties at the interface.

Table 3. Iteration profile for Test Case 3

Iterations	No. Function Eval.	% of Total	Cumulative %
0	4516	60.2	60.2
1	2028	27.0	87.2
2	566	7.6	94.8
3	362	4.8	99.6
>3	27	0.4	100

Test Case 3

For Test Case 3, the MoL scheme, with the iterative solution, was used to obtain a numerical solution of an example from Hills et al. (1989) at t = 5 days. Table 3 shows the number of function evaluations that result in 0, 1, 2, 3 or greater than 3 iterations for this test case. It is evident that the same pattern, as for Test Case 1 and 2 above, occurs with the majority of the iterative solution converging within 2 iterations. Although not explicit shown here, relative difference between the MoL iterative solution and the approximate solution from Matthews et al. (2004), over the 100 cm soil profile at t = 5 days was calculated. The calculations showed that both MoL solutions are in excellent agreement with relative error being within the order of 10^-5. Note that the difference increased around the last interface at z = 80 cm, which coincides with the wetting front reaching this boundary. This suggests that there is some difference between the approximate and iterative solution of the interface boundary condition around this wetting event but is not large enough to warrant any concern. Steady state solutions from the MoL were also comparing against the steady state analytical solution of Warrick et al. (1997). Both steady state solutions are also in excellent agreement being within 7 significant figure accuracy. Interestingly, for the Hills et al. (1989) example, the efficiency of the MoL with the iterative solution of the boundary condition was slightly better than the approximate solution from Matthews et al. (2004). The MoL solution proposed in this paper resulted in the total number of function evaluations of 7499 giving a CPU time of 6.19 mins. The MoL for Matthews et al. (2004) had 18585 function evaluations and took 7.09 mins to execute. The main difference between these to MoL schemes is that the iterative scheme essentially solves for the interface independently from the ODE integrator thereby taking the problem out of the ODE system. This results in a reduced ODE system that no longer has to explicitly deal with the node at the transition between the two soils thereby reducing the total number of function evaluations.

Conclusion

This paper has shown that the θ-form of Richards Equation can be effectively used to model unsaturated flow in layered soils. The technique incorporated an iterative solution, based on Newton’s Method, into the MoL to explicitly solve for θ at the interface between two soils. Note that the principles of the iterative solution could be used in other numerical schemes. The iterative solution was shown to be highly efficient converging mostly within 2 iterations for both fine over coarse and coarse over fine soil stratifications. The quick convergence is due to the shape of the Newton’s function (f) near the root (at f = 0) and the close proximity of the initial guess to the root. The convergence behaviour has been shown to hold as the intensity of surface flux increases and as the contrast between the two soils increase and decrease. The iterative solution has compared well against the approximate MoL solution from Matthews et al. (2004) and the steady state analytical solution of Warrick et al. (1997). Interestingly, for the case considered, the iterative MoL solution was more efficient than the approximate MoL solution completing the simulation in a short CPU time. In future work, more extensive comparisons will be made against other numerical approaches (Celia et al. 1990; Tocci et al. 1997) and the appropriateness of a primary variable switching technique to handle saturated and near saturated conditions will also be examined (Diersch and Perrochet 1999).

Acknowledgements

This work has been funded through ARC Discovery Indigenous Research Development program.

References

Celia MA, Bouloutas ET, Zarba RL (1990) A general mas-conservative numerical solution for the unsaturated flow equation. Water Resources Research 26, 1483-1496.

Colominas I, Gomez-Calvino J, Navarrina F, Casteleiro M (2002). A general numerical model for grounding analysis in layered soils. Advances in Engineering Software 33, 641-649.

Diaw EB, Lehmann F, Ackerer Ph (2001) One-dimensional simulation of solute transfer in saturated unsaturated porous media using the discontinuous finite element method. Journal of Contaminant Hydrology 51, 197-213.

Diersch H-JG, Perrochet P (1999) On the primary variable switching technique for simulating unsaturated-saturated flows. Advances in Water Resources 23, 271-301.

Heilig A, Steenhuis TS, Herbert SJ (2003) Funneled flow mechanisms in layered soil: field investigations. Journal of Hydrology 279, 210-233.

Hills RG, Porro I, Hudson DB, Wierenga PJ (1989) Modelling one-dimensional infiltration into very dry soils 1. Model development and evaluation Water Resources Research 25, 1259-1269.

Matthews CJ, Braddock RD, Sander GC (2004) Modeling flow through a one-dimensional multi-layered soil profile using the Method of Lines. Environmental Modeling and Assessment 9, 103 – 113.

Ratherfelder K, Abriola LM (1994) Mass conservative numerical solutions of the head-based Richards equation. Water Resources Research 30, 2579-2586.

Romano N, Brunone B, Santini A (1998) Numerical analysis of one-dimensional unsaturated flow in layered soils. Advances in Water Resources 21, 315-324.

Schiesser WE (1991) ‘The Numerical Method of Lines: Integration of Partial Differential Equations.’ (Academic Press Inc: San Diego).

Srivastava R, Yeh T-CJ (1991) Analytical solution for one-dimensional, transient infiltration toward the water table in homogeneous and layered soils. Water Resources Research 27, 753-762.

Tocci MD, Kelley CT, Miller CT (1997) Accurate and economical solution of the pressure-head form of Richards’ equation by the method of lines. Advances in Water Resources 20, 1-14.

van Genuchten, MTh (1980) A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal 44, 892-898.

Warrick AW, Wierenga PJ and Pan L (1997) Downward water flow through sloping layers in the vadose zone: analytical solutions for diversions. Journal of Hydrology 192, 321-337.