Simulation of Water and Contaminant Transport Through Vadose Zone - Redistribution System

Movement of water in vadose zone, mainly focusing on infiltration and percolation that involves percolation of water under gravity from soil surface and redistribution which is the capillary rise of water movement upwards, is presented. In the global hydrologic cycle, 76% of the precipitating water enters the soil via percolation-infiltration, which leads to the downward movement of water (L’vovich 1974). The water used by natural processes, can move downwards due to infiltration and lift from groundwater table during natural redistribution process. The forecasting of water movement in unsaturated infiltrationredistribution system is linked between soil hydraulic properties and hydrologic condition of natural surface water system. The understanding of water movement processes associated with infiltration and redistribution has a number of practical applications. One such application is to predict the fate and transport of materials through soil including nutrients, organic carbon and microbes under natural processes, which in turn will help in developing appropriate management plans for irrigation, fertiliser application and waste disposal on land.


Introduction
Movement of water in vadose zone, mainly focusing on infiltration and percolation that involves percolation of water under gravity from soil surface and redistribution which is the capillary rise of water movement upwards, is presented.In the global hydrologic cycle, 76% of the precipitating water enters the soil via percolation-infiltration, which leads to the downward movement of water (L'vovich 1974).The water used by natural processes, can move downwards due to infiltration and lift from groundwater table during natural redistribution process.The forecasting of water movement in unsaturated infiltrationredistribution system is linked between soil hydraulic properties and hydrologic condition of natural surface water system.The understanding of water movement processes associated with infiltration and redistribution has a number of practical applications.One such application is to predict the fate and transport of materials through soil including nutrients, organic carbon and microbes under natural processes, which in turn will help in developing appropriate management plans for irrigation, fertiliser application and waste disposal on land.

Infiltration and redistribution system
The prevailing weather condition can potentially affect the amount of water input to soil.The infiltration rate consists of (1) water input rate, which is the rate of water that arrives at soil surface due to rain, natural and artificial applications and (2) infiltration capacity, which is the maximum rate at which percolating water migrates through the soil pore.In general, water input rate responds to seasonal climatic variations and to any recharge from natural or artificial conditions.The infiltration capacity rate depends on the soil texture and soil hydraulic properties.The infiltration process can be separated into three categories (1) no ponding (2) saturation from above and (3) saturation from below.The "no ponding" refers to the condition, when water input rate is less than or equals to the infiltration capacity.The "saturation from above" presents the condition, when water input rate exceeds the 1    1) and ( 2) contains Darcy's velocity, q z , which is given by (the negative sign means downward flow) (Huyakorn et al. 1984): where z q is Darcy's velocity in vertical direction [LT -1 ].By inserting a series of tensiometers in different parts of a drainage field, the profiles of pressure head can be observed as shown in Figure 1.The negative pressure head determined in the unsaturated soil layer is due to suction head, zero pressure head occurred at the interface of groundwater table and the positive pressure head in the saturated soil layer is due to hydraulic head (gravitational head plus pore-water pressure head) (Fredlund and Rahardjo 1993).
, water moves downward with suction head (case 3).The upward flow is yielded (case 4), if Darcy's velocity is greater than zero.However, in most unsaturated zone cases, the water will experience a downward flow, these are either case 2 for evaporation or case 3 for irrigation (Warrick et al. 1991).), pressure head   and volumetric moisture content    is defined as the hydraulic properties function, which is highly nonlinear.Typically used hydraulic properties equations in this research include those of Brooks and Corey (1964), Haverkamp et al. (1977), van Genuchten (1980) and Saxton et al. (1986).Brooks & Corey (1964) hydraulic properties function was firstly derived by modifying the statistical porewater interaction models of Childs and Collis-George (1950) (cited in Fredlund & Rahadjo 1993).The Brooks & Corey (1964) model is: where a  is the air entry suction pressure head [L], S  and r  are the saturated and residual volumetric water content, respectively [unitless].M  and N ~are empirical constants.Haverkamp et al. (1977) fitted the properties of homogeneous soil in unsaturated conditions by the least square method and proposed the following empirical equations (HV equations): where A, ,  and  are the curve fitting coefficients [unitless].
van Genuchten (1980) derived the hydraulic properties equations based on the equation of Brooks & Coley (1964) and proposed the following empirical equations for hydraulic properties (VG equations): where a is the soil water retention function [L -1 ], q and p are the empirical parameters, The possible values for the coefficients presented in VG equations are given in Table 1.The coefficients are sorted by soil textures in accordance with USDA textural classes (Carsel et al. 1988).

www.intechopen.com
The hydraulic properties are estimated from the soil texture using a method generalised by Saxton et al. (1986).The textural class is assessed according to the USDA system.The water retention curve is fitted with linear regression and the formulations are in S.I. unit as follows: 1) When the hydraulic pressure is between 10 and 1500 kPa (or 1.02 to 15.3 m H 2 O), then the expression for  is given by:     The coefficient of relative permeability was obtained from the relationship between the matrix suction and volumetric moisture content.The equation was given as follows:

Numerical solution for the unsaturated flow equation
The finite discretionary scheme is given in Figure 4.The number of nodes in the system is assigned sequentially to the flow direction.As the hydraulic properties model can be applied to modify Richards' equation between the  based to  based formulations.Generally,  based Richards' equation is preferred over  based equation as it is possible to accurately measure the pressure head using tensiometers.Several previous works (Wang & Anderson 1982, Coley 1983, Segerlind 1984, Paniconi et al.1991) have assumed the approximate solution to Richards' equation as: where N j (z) is the shape function [unitless],   j t


is the unknown coefficient with corresponding to the value of nodal pressure head [L] and the subscript "j" is defined to denote a nodal sequence.
Richards' equation (Eq. 1) is now written in the form of   0 L   as follows.  The weighting function, N i was assigned the same value as j N .Applying Galerkin's finite element method yields (Bunsri et al. 2009): Substituting Eq. 13 into Eq.14, yields: The numerical solution of Richards' equation is (Bunsri et al. 2009): The algebraic matrix systems were defined as follows (Bunsri et al. 2009): where , and Using Eq. 3, the vector matrix   i E could be written in the form of Darcy's flux boundary condition as: The time derivative approximation at a particular node " j " can be explained by (Wang & Anderson 1982, Paniconi et al. 1991, Ségol 1993): where   / t   is a column matrix, which refers to time-dependent hydraulic pressure head  " can be simplified with the vector symbol "    ".The iterative scheme obtained using a single Picard's method is given by (Wang & Anderson 1982, Paniconi et al. 1991, Ségol 1993): Celie et al. ( 1990) estimated the pressure head profiles using Richards' equation.The numerical model is developed using both finite difference model (FDM) and finite element method (FEM).The model developed using FEM has given an oscillation of pressure head that is near to infiltration front.This oscillation pressure head could be reduced by applying a diagonal time matrix (mass lumping technique).The net water balance between inflow and outflow of soil pore water at any node and time step t (MB) is defined based on a diagonal time matrix.This aspect is investigated to evaluate the existence of any errors within the calculation process if existed.The water balance equation for the finite element technique is presented as follows (Celie et al. 1992, Ségol 1993): where 0 q and E q are the boundary fluxes associated with 0 z and L z , respectively [L T -1 ].
Symbols "i" and "t" count for the sequence of nodes and time steps, respectively.Subscripts "0", "E" refer to the upper boundary of downward flow and the number of elements, respectively.Superscripts "0", "t" and "t+1" refer to the initial, previous and current iteration time step, respectively.

Boundary condition for infiltration system
The initial condition was defined as the starting condition of the system at t=0.The boundary conditions were classified into upper and lower boundaries that were located on the top and bottom of the considered system (Huyakorn & Pinder 1983).The upper boundary was the condition at the discharge point and the lower boundary was the condition at the water table or the column base.The boundary and initial conditions for the solute transport model and Richards' equation included the known concentration of contaminant and pressure head, respectively (Bear 1979, Huyakorn & Pinder 1983).The change of infiltration is mainly controlled by intensity of infiltration and surface soil properties.The infiltration rates are determined by: (1) Rate of water approaching above soil surface via rainfall, snowmelt, irrigation, natural or artificial recharge and depth of ponding on the surface; (2) Fully saturated hydraulic conductivity at soil surface; (3) Degree of saturation at soil surface when the infiltration begins; (4) Inclination and roughness of topsoil; and (5) Chemical characteristics of top soil and physical and chemical properties of water.
The infiltration rate;   f t is computed as (Nassif & Wilson 1975): where i is natural precipitation or application rate and   q t is the rate of surface runoff.The natural infiltration recharge can be technically determined using Green-Ampt approach (Green & Ampt 1911).In the case of infiltration coming from runoff or rainfall recharge, the percolation rate is based on top soil properties.It can be classified into 2 categories as: 1) Water input rate is less than saturated hydraulic conductivity (i<K z ).At initial stage (t=0), the hydraulic conductivity is defined as 0 z K and it increases to K z at the specified time; t w .At this stage, water will percolate and it is stored until the soil pore is fully filled.
2) Water input rate is greater than saturated hydraulic conductivity (i>K z ).This process is observed at an early stage of infiltration in which the excess water cannot be transmitted downwards.The maximum water content and the hydraulic conductivity are limited at S  and z K , respectively.Therefore, when the soil surface reaches saturation, ponding will occur or in the case of a hilly area, overland flow will take place.However, the corresponding equation for this case cannot give a valid value of f(t) as there is no well developed relationship for ponding.The percolation of water during infiltration process with no ponding can be expressed as: where 0  is the initial soil moisture content [unitless], p f is the percolation rate, w  is the pressure head at the wetting front, I is the cumulative infiltration and S z is the wetting front distance for Green-Ampt model.Schmid (1990) has modified Green-Ampt approach with a Taylor series expansion, and obtained the explicit approximation for a function of the infiltration rate.
where f  is the effective tension at the wetting front     where f z is the wetting front distance that is located between soil surface and the bottom of the wetted soil layer.Variable  is the shape coefficient of Brook and Corey (1964) hydraulic properties model,   and  is a constant.

Boundary condition for redistribution system
The boundary conditions for a redistribution system can be evaluated according to the physical model of capillary rise.The zone of negative pressure is observed within the depth of capillary height.The capillary height is the height of the water level inside the capillary tube.The capillary height could be estimated using (Freudlune & Rahardjo): where w U is the hydraulic pressure at the capillary height h is the capillary height [L] and g is the gravitational acceleration [L T -2 ].The physical model of capillary pressure force in unsaturated soil is presented in Figure 5.If matric suction (U a -U w ) or capillary height (h C ) is plotted against pore radius (r) on a Log-Log plot, a linear relationship is expected.Further details can be found in Fredlune & Rahardjo (1993).The variables " S ", " S r "and "  " refer to the surface tension, the radius of capillary tube and the contact angle, respectively.The variable " a U " is the atmospheric pressure (guage) that is normally taken as 0 cm H 2 O.

Governing equation for solute transport
In practice, the objective is not only to predict the movement of water in the vadose zone, but also to determine the movement of reactive and non-reactive contaminants through the soil pores.The developed model can predict the fate and transport of reactive constituents such as phosphate, nitrate, organic carbon compounds and microbes.The governing equation for multi-component transportation of contaminants in porous media under variable saturation conditions could be expressed in a general form as follows (Schnoor 1996).D is an effective molecular diffusion coefficient [L 2 T -1 ],  is a coefficient relating to tortuosity [unitless] and i v is the average linear velocity in the vertical direction;

Nitrogen and organic carbon compounds
The total rate of organic carbon compounds (substrate) utilisation, S r was a combination of the substrate utilisation rate due to aerobic and nitrate respiration (Widdowson et al. 1988).

SS O
rrand SN r are the total substrate utilisation, substrate utilisation under aerobic respiration and substrate utilisation under nitrate respiration, respectively [T -1 ].Using a modified Monod's equation, the substrate utilisation rates could be derived as follows (Widdowson et al. 1988): . By referring to the rate of kinetic reaction (Eq.31a and 31b) for organic carbon compounds biodegradation, the equation for transport of organic carbon compounds is written as (Schnoor 1996): The nitrate transport equation can be formulated as follows (Schnoor 1996 (Harvey et al. 1984cited in Widdowson et al. 1988).

Phosphate phosphorus compounds
Phosphorus adsorption on soil was formulated as follows (Shah et al. 1975 C is the concentration of phosphorus in liquid phase that is in equilibrium with the concentration of phosphorus in the solid phase [M M -1 ] and * ww PP CC  is a driving force for transferring phosphorus from liquid to solid phase [M L -3 ].The Langmuir isotherm equation is the best fit for phosphorus adsorption (Shah et al. 1975).Variable * w P C is modified into Langmuir isotherm equation, yielding (Schnoor 1996, Watts 1997): where L k and M k are the coefficients (Langmuir rate constant) [unitless] and the maximum phosphate adsorption capacity on soil [L 3 M -1 ], respectively.The soluble phosphate, P C can be assumed as soil pore water, which can be presented as volumetric portion of moisture (Fetter 1992).The transport of soluble phosphate is obtained by: where  is a retardation factor

E. coli
The overall reaction rates of microbial kinetics were the summation of production, maintenance, decay, adsorption and desorption.The retardation equation was given as (Zysset et al. 1994): where n is the fraction of aqueous volume and biofilm in total volume (equals to porosity) [unitless].The concentration of E. coli relates to the substrates consumed.Substrate utilisation during metabolisation processes is defined using the first order Monod's kinetics equation as follows (Zysset et al. 1994). where   , which involves assimilation and dissimilation of microbes and t bb CC   .A general form of the mathematical model for contaminant transport coupling retardations could be written as follows (Huyakorn et al. 1985): The retardation factor,  equals 1.0 for organic carbon and nitrate compound transport equations.The biodecay factor (involving growth of microbes),  equals 0.0 for phosphate compounds transport equation.Only microbial transport contains all of these factors.

Numerical solution for solute transport equation
The approximate solution of contaminant concentration at any nodes and time "t" is defined as (Segerlind 1984, Huyakorn et al. 1985, Clement et al. 1998): where The mathematical model presented in Eq. 40 can be modified as follows: The integral form of Eq. 42 is given as follows: where   L , 0 is the extent of the vertical direction (one dimension) domain.The subscripts "i" and "j" denoted the sequence of elements in the domain as presented in Figure 4. Substituting Eq. 42 into 43, gives: The numerical solution of Eq. 44 is governed as: The algebraic matrix systems are defined as follows: The initial concentration in the entire domain   0,L at time t = 0 is defined as follows (Bear 1979, Huyakorn & Pinder 1983, Huyakorn, et al. 1985, Ségol 1993, Clement et al. 1998): where The boundary concentration on the edge of domain   0,L at time "t" is defined using the Dirichlet boundary condition (Bear 1979, Huyakorn & Pinder 1983, Huyakorn, et al. 1985, Ségol 1993, Clement et al. 1998).
The specific dispersive flux on the edge of domain   L , 0 at time "t" is employed using Neumann boundary condition.The dispersive flux was defined as z C   / , (Bear 1979, Huyakorn & Pinder 1983, Huyakorn, et al. 1985, Ségol 1993, Clement et al. 1998).
For rainfall infiltration 0 For seepage flow 0 where D C q and T C q are the portion of the boundary flux attributable to concentration due to dispersion [M L -1 T -2 ] and the portion attributable to total concentration [M L -1 T -2 ], respectively.

Model applications
The simulation results of solute transport in unsaturated infiltration-redistribution system is an active research area where a variety of attempts are being made to determine the dynamics of vadose zone in relation to water and contaminant movement and quality of soil.The model presented above is used to simulate the experimental data obtained by Paniconi et al. (1991).The input parameters for this case study are given in Table 3.The Neumann and the Dirichlet conditions were applied to the upper and lower boundary, respectively.Both the experimentally observed and simulated pressure head and moisture content profiles are presented in Figures 6 and 7, respectively.The developed models predicted the experimental observations of Paniconi et al. (1991) very well.This implied that the developed model is robust, and that it could effectively predicte the water movement in the infiltration-redistribution sytems.On the other hand, the simulation results from the solute transport model could not be compared with experimental observations due to lack of appropriate data.However, a series of simulations are carried out as shown in Figures 8 through 11.The kinetic rate constants of each contaminant are assumed based on the case study of contaminants movement in sandy soil near Perth (McArthur & Bettenay 1964, cited in Whelan & Barrow 1984).The input parameters used in model are presented in Table 4.

Domain
Column with depth of 10 m (saturation maintained at the base of the column).Table 3. Input parameters for water movement model (Paniconi et al. 1991).
Dispersivity; T  (cm) 3.25 (organic carbon), 7.54 (nitrate), 1.57 (phosphate) and 7. Figures 8, 9, 10 and 11 present the simulations of nitrate, organic carbon, phosphate and E.coli transport in infiltration-redistribution system.The simulation results reveal that the contaminants could reach the groundwater table over a longer period.Also, it appeares that the top soil can remove substantial amount of contaminants.Particularly, the organic carbon and E.coli, are removed within the few centimetres of the top soil layer.On the other hand, phosphate can move downwards to a depth of 2 metres.This might relate to phosphate adsorption capacity, which is relatively low in the sandy soil.However, to make proper assessment of contaminant transport and its potential contamination of groundwater longer periods of simulation are required.At this stage it can be concluded that the contaminant transport model presented earlier could be used to predict the contaminant transport within the unsaturated (vadose) zone.

Conclusion
A comprehensive model for predicting the movement of water and contaminants through unsaturated soil is presented.The models were developed based on Richards equation and mass balance relationships.Also, Finite Element Method (FEM) based solution techniques www.intechopen.comwere developed for obtaining numerical solution to the models developed.The water movement simulation model predicted fairly well the water movement observed through column studies.Thus, the water movement model can be used for predicting the water movement through unsaturated zone with confidence.However, it is important to use site specific input parameters for reliable results.In the case of contaminant transport model, it could not be tested with experimental observations due to lack of appropriate data.However, some of the simulation results for organic carbon, nitrate, phosphate and E.coli appear to indicate that the model is able to predict the contaminant transport through soil.It is recommended that the model is further tested with experimental and/or field data.The water can carry contaminants during percolation, however there are some chemical and biological mechanisms, which can retard the migration of contaminants.The capillary force can extract water from aquifer, which is important for redistribution system.With the infiltration-redistribution system, the pore velocity of water may be reduced.This can lead to self protection of groundwater from contamination.Further, it is apparent that the relative hydraulic conductivity K z k rw is one of the critical parameters that influence the water and contaminant transport through vadose zone.In this study the influence of relative hydraulic conductivity on the water and contaminant transport was not fully investigated.Evidently, further studies are required to fully understand the influence of relative hydraulic conductivity and hence identify critical soil types that may be readily affected by the on-site waste disposal systems.

Fig. 2 .
Fig. 2. Pressure head profiles at varying Darcy's velocities (Adapted from Warrick et al. 1991)The relationship between relative hydraulic conductivity (
soil water deficit.The soil moisture content under infiltration process can be estimated usingWang et al.  (2003)  equation as follows.

C
are the concentration of adhering microbes and free swimming microbes [M L -3 ], respectively.S C is the concentration of limiting substrate in aqueous compartment [N L -3 ], y v is the stoichiometric coefficient [M N -1 ],  is an effectiveness of biofilm [unitless]. k and d k are Monod's constant for substrate utilisation in biomass [L 3 M -1 T -1 ] and the constant of decay rate [T -1 ], respectively.S k and C k are the constant desorption (detachment) and adsorption (attachment) rate [T -1 ], respectively.The unit N is the quantity of microbes involved (cfu or MPN).
carbonaceous substrate concentration [M L -3 ] and the total microbial concentration [M L -3 ], respectively.b  and b Y are the microbial maximum specific growth rate [T -1 ] and a heterotrophic microbial yield coefficient [unitless], respectively.Sb K and m k are the substrate concentration when the rate of utilisation of half of the maximum rate under aerobic condition [M L -3 ] and a biomass maintenance rate [T -1 ], respectively.The E. coli transport equation was governed as follows (Zyset et al. 1994): where www.intechopen.comSimulation of Water and Contaminant Transport Through Vadose Zone -Redistribution System 411 0

Table 1 .
Recommended empirical coefficients for VG equations

Table 2 .
Kinetic parameters for organic carbon and nitrate retardation