Assessment of water penetration problem in unsaturated soils

Assessment of water penetration problem in unsaturated soils A. Barari, M. Omidvar, A. R. Ghotbi, and D. D. Ganji Departments of Civil and Mechanical Engineering, Babol University of Technology, Babol, Iran Faculty of Engineering, University of Golestan, Gorgan, Iran Department of Civil Engineering, Shahid Bahonar University, Kerman, Iran Received: 13 April 2009 – Accepted: 20 April 2009 – Published: 8 May 2009 Correspondence to: A. Barari (amin78404@yahoo.com) Published by Copernicus Publications on behalf of the European Geosciences Union.


Introduction
Modeling of multi-phased 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.Current equations used for describing fluid flow through porous media are based mainly on semi-empirical equations first derived by Buckingham (1907) and Richards (1931).Specifically, 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 addressing 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, 2001).This is mainly due to the dependence of hydraulic conductivity and diffusivity -two important Introduction

Conclusions References
Tables Figures

Back Close
Full 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 et al., 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., 2001;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).In geo-environmental applications, Bunsri et al. (2008) solved Richards' equation accompanied 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 Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc

Printer-friendly Version
Interactive Discussion parameters involved.In the present research, two analytical methods known as Homotopy Perturbation Method (HPM) (He, 1999a(He, , 2003;;Ganji and Sadighi, 2006;Choobbasti et al., 2008;Barari et al., 2008a;Ghotbi et al., 2008aGhotbi et al., , 2008b) ) and Variational Iteration Method (VIM) (He, 1997(He, , 1999b(He, , 2006;;Ganji et al., 2007;Ganji and Sadighi, 2007;Barari et al., 2008b) have been employed to solve the problem of one-dimensional infiltration of water in unsaturated soil governed by Richards' equation.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.

Richards' equation
The basic theories describing fluid flow through porous media were first introduced by Buckingham (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 known 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)  derive Richards' equation.Herein, one-dimensional infiltration of water in vertical direction of unsaturated soil is considered and accordingly, 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): 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): This is a fundamental equation in geotechnical engineering and is used fore modeling flow of water through unsaturated soils.For instance, the two-dimensional form of the Introducing a new term D, pore water diffusivity defined as the ration of the hydraulic conductivity and the differential water capacity, the θ-based form of Richards'equation may be obtained.D can 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 solution schemes.Brooks and Corey's model on 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 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).
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 generalized Burgers' equation is also obtained for general values of k and m (Grundy, 1983).
As seen previously, the two independent variables in Eq. ( 7) are time and depth.By applying the traveling wave technique (Wazwaz, 2005;Abdoul et al., 2008;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):

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version

Interactive Discussion
The exact solution to Eq. ( 12) can be found to be: 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 conductivity term has been selected in two independent example cases as θ 2 2 and θ 3 3 and respective n-values of one and two were associated with these conductivities.Homotopy perturbation method (HPM) (He, 1999(He, , 2003;;Abdoul et al., 2008a, b;Ganji and Sadighi, 2006;Choobbasti et al., 2008;Barari et al., 2008) and variational iteration method (VIM) (He, 1997(He, , 1999(He, , 2006;;Ganji et al., 2007;Ganji and Sadighi, 2007;Barari et al., 2008) described below have been used to solve Eq. ( 12).HPM and VIM are first introduced briefly, and are then implemented in order to solve Richards' equation as described by Eq. ( 7) and supported by Eqs. ( 10) and (11).

Basic idea of He's homotopy perturbation method (HPM)
To illustrate the basic ideas of HPM, we consider the following nonlinear differential equation: Generally speaking the operator A can be divided into a linear part L and a nonlinear part N(u).Equation ( 14) can therefore, be rewritten as: By the Homotopy technique, we construct a homotopy v(r, p) : Ω×[0, 1]→R, which satisfies: or where p∈[0, 1] is an embedding parameter, while u 0 is an initial approximation of Eq. ( 14), which satisfies the boundary conditions.Obviously, from Eqs. ( 17) and (18) we will have: 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.Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version

Interactive Discussion
According to the HPM, we can first use the embedding parameter p as a "small parameter", and assume that the solution of Eqs. ( 17) and ( 18) can be written as a power series in p: Setting p=1 yields in the approximate solution of Eq. ( 17) 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 ( 22) is convergent for most cases.However, the convergent rate depends on the nonlinear operator A(v).Moreover, He (1999) 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.
(2) 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.
5 Implementation of HPM to solve Richards' equation

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. ( 25) and rearranging the resultant equation based on powers of p-terms, one has: With the following conditions: With the effective initial approximation for v 0 from the conditions of Eq. ( 29) and solutions of Eqs.(26-28) may be written as follows: v 0 (z, t) = 0.5 − 0.5 tanh(0.25z), (30) In the same manner, the rest of components were obtained using the Maple package.

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.
A homotopy can be constructed as follows: Substituting v=v 0 +pv 1 + . . . in to Eq. ( 35) and rearranging the resultant equation based on powers of p-terms, one has: Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version

Interactive Discussion
In the same manner, the rest of components were obtained using the Maple package.
According to the HPM, we can conclude that: Therefore, substituting the values of v 0 (z, t), v 1 (z, t), v 2 (z, t) from Eqs. (40-42) in to Eq. ( 43) yields: θ(z, t) = (0.5 − 0.5 tanh(0.3333333z)) 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.Introduction

Conclusions References
Tables Figures

Back Close
Full To solve the Eq. ( 12) by means of VIM, one can construct the following correction functional: 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.25z) . (49) Using the above variational formula (48), we have: Substituting Eq. ( 49) in to Eq. ( 50) and after simplifications, we have:

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 (56) Using the above variational formula (55), we have: Substituting Eq. ( 56) in to Eq. ( 57) and after simplifications, we have: and so on.In the same way the rest of the components of the iteration formula can be obtained.Figs. 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 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.
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, nonswelling soils 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 B, f (r) and Γ are a general differential operator, a boundary operator, a known analytical function and the boundary of the domain Ω, respectively.