Analytical and Numerical Solutions of Richards' Equation with Discussions on Relative Hydraulic Conductivity

Hydraulic conductivity is of central importance in modelling both saturated and unsaturated flow in porous media. This is because it is central to Darcy's Law governing flow velocity and Richards' equation that is often used as the governing partial differential equation (PDE) for unsaturated flow. When doing numerical modelling of groundwater flow, two dominant challenges regarding hydraulic conductivity are heterogeneous media and unsaturated flow.


Introduction
Hydraulic conductivity is of central importance in modelling both saturated and unsaturated flow in porous media. This is because it is central to Darcy's Law governing flow velocity and Richards' equation that is often used as the governing partial differential equation (PDE) for unsaturated flow. When doing numerical modelling of groundwater flow, two dominant challenges regarding hydraulic conductivity are heterogeneous media and unsaturated flow.  Fig. 1 shows an example of soil layers full of heterogeneities that must be approximated in some way. Fig. 2 shows an idealization of a two-dimensional (2-D) cross section of a levee. Several layers representing different soil types are shown here. It is important to note that each layer is represented by a constant value of horizontal and vertical hydraulic conductivity rather than, for instance, a statistical variation. This is often done in numerical models and will be implemented in this work. The hydraulic conductivity values for sand and gravel are two to four times those of the silt and clay.

Heterogeneous media
www.intechopen.com An additional complexity is added in this problem by inserting the slurry wall. This type of wall is typically much less pervious than the surrounding soil, creating further stress on the computational model. This is because the numerical solution that is usually done requires a solution of a system of simultaneous, linear equations. The greater the span of orders of magnitude of hydraulic conductivity, the more challenging the solution of this system becomes.

Unsaturated flow
The last major concern and challenge discussed in this chapter regarding hydraulic conductivity with regard to computational and analytical solutions is unsaturated flow. Fig.  3 shows the location of the phreatic surface for steady-state conditions. The phreatic surface is where the ground goes from fully saturated when the soil voids are completely filled with water to partially saturated voids in the soil matrix. Above this phreatic surface, hydraulic conductivity is often modelled by where k = hydraulic conductivity of a given soil type s k = hydraulic conductivity for saturated soil r k = relative hydraulic conductivity for unsaturated soil r k is set to 1 in the saturated zone, but varies with the pressure head ( h ) in the unsaturated zone. There are many expressions for relative hydraulic conductivity in the literature and practice. Some of these will be discussed later in this chapter.

Obtaining computational results
A discretization of the flow region must be done to do the numerical analysis. Many techniques are available, but in this chapter, the finite element method (Cook, 1981) will be emphasized. Fig. 4 shows a zoom of the finite element mesh for the 2-D levee cross section given in Fig. 2 consisting of triangular elements. Define the total head as where h = pressure head φ = total head z = z coordinate or elevation Then equipotentials or total head contours can be used as a good way to visualize the data computed at each node of the mesh. Fig. 5 shows this type of plot for the levee example.
www.intechopen.com  Finally, using Darcy's Law for a homogeneous medium, where v = flow velocity www.intechopen.com A plot of velocity vectors for the levee cross section can be computed and plotted (see Fig.  6).

Relative hydraulic conductivity
One common way of representing relative hydraulic conductivity is using the van Genuchten expression (van Genuchten, 1980). First, A simpler but less useful expression for relative hydraulic conductivity is the Gardner formulation (Gardner, 1958), where α = parameter based on soil type Eq. 6 is shown here because this simpler equation is needed in the derivation of analytical solutions given later in this chapter. Regardless of the middle part of the curves, all relative hydraulic conductivity equations go from 1 at 0 h = to near 0 for negative values of h . In all these discussions, pressure head is greater than zero for saturated flow, equal to zero at the phreatic surface, and less than zero in the unsaturated zone.

Richards' equation
A common way of characterizing unsaturated flow is Richards' equation (Richards, 1931). A general version of this equation is where K = hydraulic conductivity tensor For a homogeneous, isotrophic medium, K becomes k times the identity matrix, so after using Eqs. 1 and 2, Eq. 7 becomes, Eq. 8 will be used for deriving the analytical solutions. The fact that r k is a function of h creates significant difficulty both for solving this problem numerically and deriving analytical solutions since now Eq. 8 is often severely nonlinear.

Analytical solutions
Analytical solutions are an excellent tool for checking numerical programs for accuracy. In these derivations, hydraulic conductivity plays an important role. The challenge is finding a form of relative hydraulic conductivity such that the nonlinear Richards' equation can be converted from a nonlinear to a linear form. The derivations presented here are mirrored after those presented earlier (Tracy, 2006(Tracy, , 2007 because they lend themselves to onedimensional (1-D), 2-D, and three-dimensional (3-D) solutions. First, 1-D and 2-D analytical solutions will be derived, and then numerical finite element solutions highlighting accuracy for different representations of relative hydraulic conductivity will be investigated. Fig. 7 shows the 1-D problem that will be considered in detail. A column of soil of height, L , is initially dry until water begins to infiltrate the soil. A pool of water at the ground surface is then maintained holding the pressure head to zero. This is known as the 1-D Green-Ampt problem (Green & Ampt, 1911). Fig. 7. A view of a 1-D column of soil that is initially dry until water is applied at the top of the ground surface from rainfall.

www.intechopen.com
This problem is challenging numerically because the change in relative hydraulic conductivity is so dramatic, as it goes from small to one. There are several steps that are involved in the derivation for this problem, and they will now be summarized. 1. Provide a function of relative hydraulic conductivity and moisture content as a function of pressure head. 2. Establish initial and boundary conditions. 3. Perform a change of variables to linearize Richards' equation. 4. Solve this new PDE for the steady-state solution. 5. Obtain yet another PDE using a second change of variables. 6. Use separation of variables. 7. Use Fourier series to solve the current PDE. 8. Transform back to the original variables.

Relative hydraulic conductivity and moisture content
Gardner's equation (Eq. 6) is used for relative hydraulic conductivity, and moisture content is given by where d θ = moisture content when the soil is dry s θ = moisture content when the soil is saturated Rather than use the van Genuchten expression for e S , a simpler version is used (Warrick, 2003) as follows: This equation is more limiting in actual practical application, but it allows easier derivation of the analytical solution. It is certainly good enough to test different computational strategies in computer programs.

Initial and boundary conditions
The initial conditions are that the soil is dry. Thus, where d h = the pressure head when the soil is dry www.intechopen.com

Change of variables
The 1-D version of Eq. 8 is with initial conditions, ( ) and boundary conditions for 0 t > from Eq. 12,

Steady-state solution
The steady-state version of Eq. 19 will now be solved. It is important to note that this steadystate version now becomes an ordinary differential equation (ODE) as follows: www.intechopen.com

Separation of variables
Eq. 28 can be solved using separation of variables. ĥ will be cast into the form, is a function only of t . Substituting Eq. 30 into Eq. 28 and dividing by ΖΤ gives The only nontrivial solution occurs when the left-and right-hand sides of Eq. 31 are set to the same arbitrary constant, η . Thus, This leads to the characteristic equations, The solution for ĥ then becomes

Transform back
The last remaining task is to convert back to the original coordinates using Eqs. 14, 25, and 26. Therefore, www.intechopen.com Rainfall L a Fig. 8. A view of a 2-D cross section of soil that is initially dry until water is applied at the top

Analytical solution of a 2-D infiltration problem
The great thing about the above derivations is that they can be extended to two and three dimensions. Fig. 8 shows a 2-D cross section of a region of soil of dimensions, a × L , where a 2-D Green-Ampt problem is presented. The soil is initially dry until water is supplied such that a specified pressure head is applied at the top with pressure head set to zero in the middle and tapering rapidly to d h at 0 x = and xa = . Fig. 9 shows the function selected to achieve this for d h = -20 m, and a = 50 m.
The equation for h is now

Transient solution for ĥ
The equation for ĥ is now with initial and boundary conditions, as before, Also, as before, transforming back to the original h , A 3-D solution is done in a similar manner.

Numerical models
Hydraulic conductivity has an important role in numerical models. Many soil layers can be modelled by specifying hydraulic conductivity for the different layers. Because Richards' equation is nonlinear, the manner in which numerical models compute relative hydraulic conductivity is also important for both accuracy of the solution and the ability of the numerical algorithms to converge. When doing a 3-D Green-Ampt problem containing thousands of 3-D finite elements on a parallel high performance computing platform, the solution would not converge because of how relative hydraulic conductivity was computed inside each finite element. When the pressure head was averaged from the four nodes of each tetrahedral element and then used to compute a constant value for the relative hydraulic conductivity inside the element, the solution diverged. However, if relative hydraulic conductivity was considered to vary linearly inside each element, the solution converged quite well. Testing these different algorithms is greatly enhanced by the analytical solutions presented above. Some tests using the analytical solutions will now be illustrated.

1-D solution of the Green-Ampt problem
The 1-D version of Eq. 7 for a homogeneous, isotropic soil is A finite element/finite difference/finite volume discretization of this equation (see Fig. 10 where j = node number r k − = the relative hydraulic conductivity for the element between nodes j and 1 j − r k + = the relative hydraulic conductivity for the element between nodes j and 1 j + t Δ = time-step size n = time-step number The two ways of computing relative hydraulic conductivity inside each element will now be discussed.

2-D solution of the Green-Ampt problem
The 2-D version of Eq. 7 was solved for the problem given in Section 4.2 with the values of the parameters being the same as for the 1-D problem presented above but with the addition of a = 50 m. The model used for this computation was a transient version of Seep2D (Tracy, 1983, & Seep2D, 2011. A steady-state version of Seep2D is currently incorporated into the www.intechopen.com Groundwater Modeling System (GMS) (Jones, 1999, & GMS, 2011. The transient version is not yet available. Fig. 11 gives a color contour plot of the error for the linearly varying relative hydraulic conductivity option for α = 0.1 m -1 for the upper, right-hand region of 10 m × 25 m. Clearly, the results match well with the analytical solution.

Summary
This chapter has shown that hydraulic conductivity plays an important role in both deriving analytical solutions and doing numerical computations. Analytical solutions for both the 1-D and 2-D Green-Ampt problem were derived and computed numerically with the results compared. The derivations are presented in such detail that others can do additional solutions as well. Varying relative hydraulic conductivity linearly within each finite element not only makes the nonlinear convergence algorithm more robust, but it also produces more accurate answers than when it is considered constant inside each finite element.

Acknowledgment
This work was supported in part by a grant of computer time from the DoD High Performance Computing Modernization Program. There are several books on broad aspects of hydrogeology, groundwater hydrology and geohydrology, which do not discuss in detail on the intrigues of hydraulic conductivity elaborately. However, this book on Hydraulic Conductivity presents comprehensive reviews of new measurements and numerical techniques for estimating hydraulic conductivity. This is achieved by the chapters written by various experts in this field of research into a number of clustered themes covering different aspects of hydraulic conductivity. The sections in the book are: Hydraulic conductivity and its importance, Hydraulic conductivity and plant systems, Determination by mathematical and laboratory methods, Determination by field techniques and Modelling and hydraulic conductivity. Each of these sections of the book includes chapters highlighting the salient aspects and most of these chapters explain the facts with the help of some case studies. Thus this book has a good mix of chapters dealing with various and vital aspects of hydraulic conductivity from various authors of different countries.