Numerical analysis of Richards ’ problem for water penetration in unsaturated soils

Unsaturated flow of soils in unsaturated soils is an important problem in geotechnical and geo-environmental engineering. Richards’ equation is often used to model this phenomenon in porous media. Obtaining appropriate solution to this equation therefore provides better means to studying the infiltration into unsaturated soils. Available 5 methods for the solution of Richards’ equation mostly fall in the category of numerical methods, often having restrictions for practical cases. In this research, two analytical methods known as Homotopy Perturbation Method (HPM) and Variational Iteration Method (VIM) have been successfully utilized for solving Richards’ equation. Results obtained from the two methods mentioned show a remarkably high precision in the 10 obtained solution, compared with the existing exact solutions available.


Introduction
Modeling water flow through porous media presents an important problem of practical interest for geotechnical and geo-environmental engineering, as well as many other areas of science and engineering.Study of this phenomenon requires proper formulation of the governing equations and constitutive relations involved.Currently, equations used for describing fluid flow through porous media are based mainly on semi-empirical equations first derived by Buckingham (1907) and Richards (1931).Despite limitations and drawbacks, Richards' equation is still the most widely used equation for modeling unsaturated flow of water through soil (porous media) (Hoffmann, 2003).Due to the importance and wide applications of the problem, many researches have been devoted in the past to proper assessment of different forms of Richards' equation.Both analytical and numerical solutions have been investigated in the literature.Analytical solutions to Richards' equation are rather scarce and are generally limited to only special cases (Ju and Kung, 1997;Arampatzis et al., 2001).This is mainly due to the dependence of hydraulic conductivity and diffusivity -two important parameters in the equation -on moisture content, combined with the non-trivial forcing conditions that are often encountered in engineering practice (Ju and Kung, 1997;Arampatzis et al., 2001;Kavetski et al., 2002).As a result, application of many numerical methods to the solution of Richards' equation with various engineering applications has been investigated in the literature.Finite element and finite difference methods have been adopted by several researchers (Clement et al., 1994;Baca et al., 1997;Bergamaschi and Putti, 1999;Milly, 1985).Mass lumping was employed in these studies to improve stability.Time stepping schemes such as the Douglas-Jones-predictor-corrector method, Runge-Kutta method and backward difference formulae should also be mentioned in this context (Kavetski et al., 2001a;Miller et al., 2005).Tabuada et al. (1995) used an implicit method and presented equations governing two-dimensional irrigation of water into unsaturated soil based on Richards' equation.The Gauss-Seidel method was then effectively used to solve the resulting equations.Ross (2003) introduced an efficient non-iterative solution for Richards' equation using soil property descriptions as proposed by Brooks and Corey (1964).In his method, Ross used a space and time discretization scheme in order to derive a tridiagonal set of linear equations which were then solved non-iteratively.Varado et al. (2006) later conducted a thorough assessment of the method proposed by Ross and concluded that the model provides robust and accurate solutions as compared with available analytical solutions (Basha, 1999).Several other iterative solutions have also been cited in the literature.One such study is that of Farthing et al. (2003) which used the well-known pseudo-transient continuation approach to solve the nonlinear transient water infiltration problem, as well as the steady-state response as governed Richards' equation.Other commonly used iterative schemes include the Picard iteration scheme (Chounet et al., 1999;Forsyth et al., 1995), the Newton and inexact Newton schemes (Jones and Woodward, 2001;Kavetski et al., 2001b;Kees and Miller, 2002) and hybrid Newton-Picard methods.Huang et al. (1996) considered the modified Picard iteration schemes and presented several convergence criteria as to evaluate the efficiency of the various iterative methods.In geo-environmental applications, Bunsri et al. (2008) solved Richards' equation accom-6361 panied by advective-dispersive solute transport equations by the Galerkin technique.Witelski (1997) used perturbation methods to study the interaction of wetting fronts with impervious boundaries in layered soils governed by Richards' equation.Through comparison with numerical solutions, Witelski concluded that perturbation methods are able to yield highly accurate solutions to Richards' equation (Witelski, 1997).
Each method mentioned above encounters Richards' equation in a certain way.Often, assumptions are made and empirical models are implemented in order to overcome difficulties in solving the equation due to high interdependence of some of the parameters involved.Analytical solutions often fit under classical perturbation methods (Kevorkian and Cole, 1996;Nayfeh, 1973;Nayfeh and Mook, 1979).However, as with other analytical techniques, certain limitations restrict the wide application of perturbation methods, most important of which is the dependence of these methods on the existence of a small parameter in the equation.Disappointingly, the majority of nonlinear problems have no small parameter at all.Even in cases where a small parameter does exist, the determination of such a parameter doesn't seem to follow any strict rule, and is rather problem-specific.Furthermore, the approximate solutions solved by the perturbation methods are valid, in most cases, only for the small values of the parameters.It is obvious that all these limitations come from the small parameter assumption.In the present study, two powerful analytical methods Homotopy Perturbation Method (HPM) (He, 2003(He, , 1999a(He, , 2006a(He, , 2000;;Barari et al., 2008a,b;Ghotbi et al., 2008a,b) and Variational Iteration Method (VIM) (He, 1997(He, , 1999b(He, , 2006b;;Sweilam and Khader, 2007;Momani and Abuasad, 2006;Barari et al., 2008c) have been employed to solve the problem of one-dimensional infiltration of water in unsaturated soil governed by Richards' equation.HPM is actually a coupling of the perturbation method and Homotopy Method, which has eliminated limitations of the traditional perturbation methods.HPM requires no small parameters in the equations and can readily eliminate the limitations of the traditional perturbation techniques.He (He, 1997(He, , 1999b(He, , 2006b) ) proposed a variational iteration method (VIM) based on the use of restricted variations and correction functionals which has found a wide application for the so-lution of nonlinear ordinary and partial differential equations.This method does not require the presence of small parameters in the differential equation, and provides the solution (or an approximation to it) as a sequence of iterates.The method does not require that the nonlinearities be differentiable with respect to the dependent variable and its derivatives.
In the next sections, Richards' equation and the relative models involved are introduced, followed by a thorough explanation of the analytical methods used to solve the equation.Illustrative examples are also given in order to demonstrate the effectiveness of the method in solving Richards' equation.The examples considered, despite involving relatively simple one-dimensional cases, represent a nonlinear problem.Despite the major advances of science in the field of solving differential equations, it is still very difficult to solve such nonlinear problems either numerically or analytically.He and Lee (2009) attributes this shortcoming to the fact that various discredited methods and numerical simulations apply iteration techniques to find numerical solutions of nonlinear problems, and since nearly all iterative methods are sensitive to initial solutions, it is very difficult to obtain converged results in cases of strong nonlinearity.

Richards' equation
The basic theories describing fluid flow through porous media were first introduced by Buckingham (1907) who realized that water flow in unsaturated soil is highly dependent on water content.Buckingham introduced the concept of "conductivity", dependent on water content, which is today known as unsaturated hydraulic conductivity (after Rolston, 2007).This equation is usually referred to as Buckingham law (Narasimhan, 2005).Buckingham also went on to define moisture diffusivity which is the product of the unsaturated hydraulic conductivity and the slope of the soil-water characteristic curve.Nearly two decades later, Richards (1931) applied the continuity equation to Buckingham's law -which itself is an extension of Darcy's law -and obtained a general partial differential equation describing water flow in unsaturated, non-swelling soils 6363 with the matric potential as the single dependent variable (Philip, 1974).There are generally three main forms of Richards' equation present in the literature namely the mixed formulation, the h-based formulation and the θ-based formulation, where h is the weight-based pressure potential and θ is the volumetric water content.
Since Richards' equation is a general combination of Darcy's law and the continuity equation as previously mentioned, the two relations must first be written in order to derive Richards' equation.Herein, one-dimensional infiltration of water in vertical direction of unsaturated soil is considered, for which Darcy's law and the continuity equation are given by Eqs. ( 1) and ( 2) respectively: and where K is hydraulic conductivity, H is head equivalent of hydraulic potential, q is flux density and t is time.The mixed form of Richards' equation is obtained by substituting Eq. ( 1) in Eq. ( 2): (3) Equation ( 3) has two independent variables: the soil water content, θ, and pore water pressure head, h.Obtaining solutions to this equation therefore requires constitutive relations to describe the interdependence among pressure, saturation and hydraulic conductivity.However, it is possible to eliminate either θ or h by adopting the concept of differential water capacity, defined as the derivative of the soil water retention curve: The h-based formulation of Richards' equation is thus obtained by replacing Eq. ( 4) in Eq. ( 3): (5) This is a fundamental equation in geotechnical and geo-environmental engineering and is used for modeling flow of water through unsaturated soils.For instance, the twodimensional form of the equation can be used to model seepage in the unsaturated zone above water table in an earth dam.
Introducing a new term D, pore water diffusivity, defined as the ratio of the hydraulic conductivity to the differential water capacity, the θ-based form of Richards' equation may be obtained.D can therefor be written as: It should be noted that both D and K are highly dependent on water content.Combining Eq. ( 6) with Eq. ( 3) gives Richards' equation as: In order to solve Eq. ( 7), one must first properly address the task of estimating D and K , both of which are dependent on water content.Several models have been suggested for determining these parameters.The Van Genuchten model (Van Genuchten, 1980) and Brooks and Corey's model (Brooks and Corey, 1964;Corey, 1994) are the more commonly used models.The Van Genuchten model uses mathematical relations to relate soil water pressure head with water content and unsaturated hydraulic conductivity, through a concept called "relative saturation rate".This model matches experimental data but its functional form is rather complicated and it is therefore difficult to implement it in most analytical solution schemes.Brooks and Corey's model on 6365 the other hand has a more precise definition and is therefore adopted in the present research.This model uses the following relations to define hydraulic conductivity and water diffusivity: Where D(θ) and K (θ) represent diffusivity and conductivity, respectively, K s is saturated conductivity, θ r is residual water content, θ s is saturated water content and α and λ are experimentally determined parameters.Brooks and Corey determined λ as pore-size distribution index (Brooks and Corey, 1964).A soil with uniform pore-size possesses a large λ while a soil with varying pore-size has small λ value.Theoretically, the former can reach infinity and the latter can tend towards zero.Further manipulation of Brooks and Corey's model yields the following equations (Witelski, 1997;Corey, 1986;Witelski, 2005): where K 0 , D 0 and k are constants representing soil properties such as pore-size distribution, particle size, etc.In this representation of D and K , θ is scaled between 0 and 1 and diffusivity is normalized so that for all values of m, ∫ D(θ)dθ=1 (after Nasseri et al., 2008).Equation ( 11) suggests that conductivity may have linear, parabolic, cubic, etc. variation with water content, associated with k values of 1, 2, 3, etc., respectively.
Several analytical and numerical solutions to Richards' equation exist based on Brooks and Corey's representation of D and K .Replacing n=0 and k=2 in Eqs. ( 10) and ( 11) yields the classic Burgers' equation extensively studied by many researchers (Basha, 2002;Broadbridge and Rogers, 1990;Whitman, 1974).The gen-eralized Burgers' equation is also obtained for general values of k and n (Whitman, 1974).
As seen previously, the two independent variables in Eq. ( 7) are time and depth.By applying the traveling wave technique (Wazwaz, 2005;He, 1997;Elwakil et al., 2004), instead of time and depth, a new variable which is a linear combination of them is found.Tangent-hyperbolic function is commonly applied to solve these transform equations (Soliman, 2006;Abdou and Soliman, 2006).Therefore the general form of Burgers' equation in order of (n,1) is obtained as (Wazwaz, 2005): The exact solution to Eq. ( 12) can be found to be: A 2 = γα 1 + n γ is an arbitrary coefficient which is selected as 1 here, following Nasseri et al. (2008).In this study, nonlinear infiltration of water in unsaturated soil has been studied using Richards' equation and by employing Brooks and Corey's model to represent hydraulic conductivity and diffusivity.The problems discussed herein may be generalized to any combination of (k, n), having different physical interpretations.However, as Witelski (1997) mentioned, the case of k = n + 1 provides certain analytical simplifications, and will therefore be considered here.The conductivity function has been selected in two independent example cases as θ 2 /2 and θ 3 /3 which corresponds to parabolic and cubic variation of conductivity with k=2 and k=3, respectively.The linear variation of conductivity, i.e., k=1 represents a linear problem and has already been thoroughly studied by Witelski (1997).The respective n values of one and two are associated 6367 with the above chosen conductivities.The values considered represent fine-textured to medium-textured soils.The resulting form of Richards' equation will therefore be: Witelski (1997) reasoned that nonlinear variations of conductivity means that steadyprofile traveling wave solutions will exist in the solutions that move with constant speed and balance conductive transport and dispersal.The profiles do not vary with time, unlike the linear problem which has a diffusive nature.

Basic idea of He's homotopy perturbation method (HPM)
Simply put, solving a nonlinear problem by means of HPM begins with a trial-function of the same order of the original problem with some unknown parameters, followed by constructing a linear differential equation whose solution is the chosen trial-function.The next step is to construct such a homotopy that when the homotopy parameter p=0, it becomes the above constructed linear equation; and when p=1, it turns out to be the original nonlinear equation.The changing process of p from zero to unity is just that of the trial-function (initial solution) to the exact solution.To approximately solve the problem, the solution is expanded into a series of p, just like that of the classical perturbation method.Generally, one iteration is enough, but the solution is always obtained with three steps at most.
To illustrate the basic ideas of HPM, we consider the following nonlinear differential equation: with the boundary conditions of where A, B, f (r) and Γ are a general differential operator, a boundary operator, a known analytical function and the boundary of the domain Ω, respectively.Generally speaking the operator A can be divided into a linear part L and a nonlinear part N(u).Equation ( 15) can therefore, be rewritten as: By the Homotopy technique, we construct a homotopy v(r, p):Ω×[0, 1]→R, which satisfies: where p∈[0,1] is an embedding parameter, while u 0 is an initial approximation of Eq. ( 15), which satisfies the boundary conditions.Obviously, from Eqs. ( 18) and ( 19) we will have: 6369 The changing process of p from zero to unity is just that of v(r, p) from u 0 (r) to u(r).In topology, this is called deformation, while L(v)−L(u 0 ) and A(v)−f (r) are called homotopy.
According to the HPM, we can first use the embedding parameter p as a "small parameter", and assume that the solution of Eqs. ( 18) and ( 19) can be written as a power series in p: Setting p=1 yields in the approximate solution of Eq. ( 18) to: The combination of the perturbation method and the homotopy method is called the HPM, which eliminates the drawbacks of the traditional perturbation methods while keeping all its advantage.The series ( 23) is convergent for most cases.However, the convergent rate depends on the nonlinear operator A(v).Moreover, He (1999a) made the following suggestions: 1.The second derivative of N(v) with respect to v must be small because the parameter may be relatively large, i.e. p→1.

The norm of L −1 ∂N
∂v must be smaller than one so that the series converges.

Basic idea of variational iteration method (VIM)
To clarify the basic ideas of VIM, we consider the following differential equation: where L is a linear operator, N is a nonlinear operator and g(t) is an inhomogeneous term.
According to VIM, we can write down a correction functional as follows: Where λ is a general lagrangian multiplier which can be identified optimally via the variational theory.The subscript n indicates the n-th approximation and u n is considered as a restricted variation, i.e. δ ũn =0.

Case 1: if n=1
To solve Richard equation when n=1 by means of HPM, we consider the following process after separating the linear and nonlinear parts of the equation.
A homotopy can be constructed as follows: Substituting v=v 0 +pv 1 +... in to Eq. ( 26) and rearranging the resultant equation based on powers of p-terms, one has: With the following conditions: v 0 (z, 0) = 0.5 − 0.5 tanh(0.25z) With the effective initial approximation for v 0 from the conditions (30) and solutions of Eqs.(27-29) may be written as follows: In the same manner, the rest of components were obtained using the Maple package.
According to the HPM, we can conclude that:

Case 2: if n=2
To solve Richard equation if n=2 by means of HPM, we consider the following process after separating the linear and nonlinear parts of the equation.
The obtained solution has been drawn in Figs. 1 and 2 for n=1 and n=2, respectively, along with the results of VIM from the following section, as well as the exact solution available.It is noteworthy from the figures that the plotted values of θ(z, t) include negative z-values as well as positive values, whereas Richards' equation usually has significant physical meaning for positive z-values.However, since there are some solutions given for infinite conditions in Richard's equation in the literature, i.e., for −∞<z<+∞ (especially for linear variations of k), solutions for negative values of z have also been presented in the figures.
6 Implementation of VIM to solve Richards' equation

Case 2: if n=2
To solve the Richards' equation by means of VIM, once again the following correction functional may be constructed: Its stationary conditions can be obtained as follows: We obtain the lagrangian multiplier: As a result, we obtain the following iteration formula: Now we start with an arbitrary initial approximation that satisfies the initial condition: θ 0 (z, t) = (0.5 − 0.5 tanh(0.3333z)) 0.5 .(57) Using the above variational formula (56), we have: Substituting Eq. ( 57) in to Eq. ( 58) and after simplifications, we have: 6377 and so on.In the same way the rest of the components of the iteration formula can be obtained.Figures 1 and 2 show the results obtained from VIM and HPM along with the exact solution, revealing a high level of agreement between the three results shown.It should be noted that Richards' equation has significant physical meaning for positive z-values.However, since there are some solutions given for infinite conditions in Richard's equation in the literature, i.e., for −∞<z<+∞ (especially for linear variations of k), solutions for negative values of z have also been presented in the figures.It is evident from the curves plotted that the exact solution and the obtained solutions from HPM and VIM almost completely overlay each other and the level of agreement between the results is therefore excellent.

Concluding remarks
In this study, Analytical solution to Richards' equation was explored using Homotopy Perturbation Method (HPM) and Variational Iteration Method (VIM).Richards' equation is used for modeling infiltration in unsaturated soils.HPM and VIM have been successfully utilized for solving Richards' equation.Illustrative examples proved the high accuracy of the results obtained using HPM and VIM.Compared to the methods presented in the literature, the homotopy perturbation method is an asymptotic method in the sense that no convergence proof is needed for the solutions obtained, i.e. a converged solution is achieved within three iterations at the most, which means the method is quite robust in achieving the solution.Although the examples considered for Richards' equation herein involved one-dimensional cases, the success of the methods used are still appraisable, as it should be noted that the equation is still a nonlinear one even in its one-dimensional form.Following the success observed in one-dimensional cases, it is now hoped that more complex situations may also be addressed in the future.It can be concluded that HPM and VIM may be effectively used for solving Richards' equation which covers some important class of problems in geotechnical and geo-environmental engineering.
and so on.In the same way the rest of the components of the iteration formula can be obtained.Figs.
(1) and ( 2   6384 Compared to the methods presented in the literature, the homotopy perturbation method is an asymptotic method in the sense that no convergence proof is needed for the solutions obtained, i.e. a converged solution is achieved within three iterations at the most, which means the method is quite robust in achieving the solution.Although the examples considered for Richards' equation herein involved one-dimensional cases, the success of the methods used are still appraisable, as it should be noted that the equation is still a nonlinear one even in its one-dimensional form.Following the success observed in one-dimensional cases, it is now hoped that more complex situations may also be addressed in the future.It can be concluded that HPM and VIM may be effectively used for solving Richards' equation which covers some important class of problems in geotechnical and geoenvironmental engineering. ) show the results obtained from VIM and HPM along with the exact solution, revealing a high level of agreement between the three results shown.It should be noted that Richards' equation has significant physical meaning for positive z-values.However, since there are some solutions given for infinite conditions in Richard's equation in the literature, i.e., for z −∞ < < +∞ (especially for linear variations of k), solutions for negative values of z have also been presented in the figures.It is evident from the curves plotted that the exact solution and the obtained solutions from HPM and VIM almost completely overlay each other and the level of agreement between the results is therefore excellent.

Figure 1 :
Figure 1: Plot of θ(z) for different values of time (t=0,1,3,5) considering n=1 -solid lines represent results from VIM and HPM, while dashed lines depict the exact solution.

Figure 2 :
Figure 2: Plot of θ(z) for different values of time (t=0,1,3,5) considering n=2 -solid lines represent results from VIM and HPM, while dashed lines depict the exact solution.