2. LITERATURE REVIEW

A review of the literature revealed that virtually no work has been published on a combined heterogeneous vadose zone transport model and saturated zone mixing model. Some research has been conducted to estimate the impact due to the mobilization and migration of contaminant in the vadose zone on the underlying groundwater resource (Ravi and Johnson 1993; Lee, 1996). Past research work focused mainly on homogeneous vadose zone leaching impact from agricultural applications (Kool and van Genuchten, 1991; Nofziger et al., 1988; Nofziger et al., 1994; Stevens et al., 1989).

2.1 Vadose Zone

Organic contaminants are transported or changed via several possible pathways including advection, dispersion, sorption, decay of the contaminant, volatilization, plant uptake, degradation and/or transformation reactions, and retention in the soil as bound residues (Jury et al., 1983; Jury, 1986; Lee, 1996).

Conceptualization of all of the above possible processes that determine the fate and transport of organic contaminants is the first step in the development of a model; however no model can simulate these processes without simplifying assumptions and a necessary digression from mechanistic, deterministic, or physically descriptive approaches to functional, stochastic, or empirical approaches to model development.

2.1.1 Structure of Vadose Zone

The early literature recognized three divisions within the vadose zone: the capillary fringe, the intermediate belt (gravitational zone) and the belt of soil water (root zone) (Meinzer, 1923). According to Meinzer, the capillary fringe is "a zone in which the pressure is less than atmospheric, overlying the zone of saturation and containing capillary interstices some or all of which are filled with water that is continuous with the water in the zone of saturation but is held above that zone by capillarity acting against gravity." That is, the capillary fringe is a saturated zone above the water table where water is affected by capillary forces.

The uppermost belt, or belt of soil water, is "that part of the lithosphere immediately below the surface, from which water is discharged into the atmosphere in perceptible quantities by the action of plants or by soil evaporation." This definition recognizes that plants, for the most part, extract water from a portion of the soil, the "root zone," near the surface.

The intermediate belt is "that part of a zone of aeration (the vadose zone) that lies between the belt of soil water and the capillary fringe." The intermediate belt is distinguished mainly by the fact that something must be between the root zone and the capillary fringe (Figure 2.1).

As shown in Figure 2.1, volumetric moisture content, q, is defined as the volume of water per bulk volume of soil sample. After a rather rapid decrease from saturation, the moisture content in the intermediate belt may remain fairly constant. Field capacity is a term used to represent this "constant" moisture content. The moisture content in the soil water belt decreases rapidly from the field capacity due to the extraction of water by plant roots and to direct evaporation at the soil surface.

Figure 2.1 Schematic diagram of capillary fringe interface between vadose (unsaturated) zone and saturated zone. The distribution of moisture in the vadose zone and the classification of waters according to Meinzer (1923).

As shown in Figure 2.2, The movement of an organic contaminant in the vadose zone can be described within and among three different phases: (1) as a solute dissolved in water, (2) as a gas in the vapor phase, and (3) as an adsorbed compound in the solid phase (Ravi and Johnson, 1993).

Figure 2.2 Schematic diagram illustrating the contaminant pathways in the vadose zone simulated by Vadose Zone Leaching Sub-model (from Ravi and Johnson, 1993).

2.1.2 Vadose Zone Transport

Leaching

Leaching refers to the downward movement of organic contaminants dissolved into the soil solution, with the primary mechanism of chemical displacement being the redistribution of water within the soil profile (Wagenet and Rao, 1985).

Processes involved in solute transport through soil are mass flow (convection), molecular diffusion, and hydrodynamic dispersion. Carefully kept daily records of water inputs and potential evapotranspiration are thus utilize some form of the convective-dispersive transport equation, which has been derived in detail (Kirkham and Powers, 1972). The general analytical form for water flow is

where C is the solute concentration [mg/L]. D is the hydrodynamic dispersion coefficient [ft²/yr]. q is the pore water velocity [ft/yr]. x is the distance [ft] and t is time.

Gish and Coffman (1987) successfully combined the general model with a stochastic model for the movement of bromide through an established no-till cornfield. This and a series of related literature indicates that most of the variability associated with chemical dispersion can be attributed to the vertical water velocity regardless of the presence or absence of plants (Addiscott and Wagenet, 1985). The value of stochastic elements in a modeling approach to account for spatial variability of soil properties is appreciated here. However, application of the convective-dispersive equation with additional terms to describe physical, chemical, and biological interactions of organic pesticides with the soil-water regime is still limited to laboratory experiments due to data-intensive requirements of such a model (Wagenet and Rao, 1985; Ravi and William, 1998; Williams et al., 1998).

Simplified approaches to water and solute transport used in most field-scale or management-oriented models are based upon the capacity, of the soil to hold water, and "piston displacement" (Wagenet and Rao, 1985). Field capacity, saturated soil moisture content, the air-dry water content, and the water content at the wilting point are often the parameters utilized to produce water fluxes or the depth of solute leaching. It seems that fluxes, not mass or concentrations, only are feasible (Rao and Wagenet, 1984). A simple model for field-scale simulations has advantages in that it is easy to design, and necessary input can be provided from existing data bases (Rao and Jessup, 1982).

Sorption

Sorption is the degree to which chemicals adhere to soil particles. Highly sorbed chemicals to soil organic matter and clay particles tend to remain near the soil surface unless preferential flow paths are available for the transport of these soil particles. Phosphorus is highly sorbed to clay particles, so it does not leach through clayey soils. Moderately and poorly sorbed chemicals have a higher potential to reach groundwater depending on the amount and timing of water applications (Jury, 1986).

The linear and Freundlich isotherms are the main models used for pesticide adsorption-desorption on soils and sediment (Rao and Jessup, 1982). These equations are:

where S is the adsorbed-phase concentration [mg/L]. C is the solution phase concentration [mg/L]. K_d and K_f are sorption coefficients. N is an empirical adsorption isotherm constant (N<1). Where isotherms depend on rate-limited processes, the above equilibrium equations become partial differential equations (Rao and Davidson, 1980). Nevertheless, equilibrium linear isotherms may be adequate over the range of solution concentrations associated with agricultural applications of pesticides (Rao and Davidson, 1980).

The single best indicator of the sorption coefficients, K_d and K_f, is the soil organic carbon content (Wagenet and Rao, 1985 ). The large variability in K_d and K_f due to adsorption by different soils can be greatly reduced by normalizing adsorption based on the organic carbon content of the soil, K_oc (Jury, 1983).

Spatial and Temporal Variability

Spatial variability refers to changes in a soil property, such as bulk density, water content, total porosity, macroporosity, hydraulic conductivity, organic carbon content, sorption coefficients, degradation rates, etc. over the field area at a given depth, and also with depth at a given location (Wagenet and Rao, 1985). Temporal variability, or changes in a soil property at a given location with time, must also be considered in the evaluation of chemical fate in the unsaturated zone, but there is little published data on temporal variability (Wagenet and Rao, 1985).

Spatial variability is usually taken as the combined results of extrinsic and intrinsic factors (Wagenet and Rao, 1985). Extrinsic factors that contribute to variability include a number of field management treatments, such as tillage, irrigation, and pesticide and fertilizer applications. In a study of a loamy sand following tillage, bulk density and the cone index become significantly greater, and total porosity and macroporosity became significantly smaller in temporal data by position interaction for the 0 to 14 cm depth (Cassel, 1983).

Intrinsic variability refers to natural variation in physical, chemical, and biological properties of field soils (Wagenet and Rao, 1985). An example would be heterogeneous soil texture due to erosion or deposition processes.

In modeling chemical fate and movement, it is not feasible to treat three-dimensional spatial variability in a deterministic manner; instead, a stochastic approach was used to quantify spatially variable field data: Fisher statistics and geostatistics (regionalized variable theory), with the emphasis on geostatistics. It is recommended that sampling programs be designed so that geostatistics and model validation can be performed.

2.1.3   Numerical Leaching Models

VLEACH (Ravi and Johnson, 1993) and VLEACHSM 1.0a (Lee, 1996) were developed as one-dimensional vertical homogeneous vadose zone leaching models. The VLEACH simulates vertical homogeneous transport by advection in the liquid phase and by gaseous diffusion in the vapor phase. The VLEACHSM 1.0a added vertical homogeneous transport to simulate the liquid-phase dispersion.

Mathematical Discussion of VLEACH

The VLEACH treats spatial and temporal variation of contaminant concentration in solid, liquid, and gas phases. These variables are denoted as C_s(z,t), C_w(z,t), and C_a(z,t), respectively, where z and t stand for the space and time variables. The total contaminant mass initially (prior to infiltration of water) present in the soil is assumed to be dissolved into the liquid phase. This yields the following relation:

where C_w(z,0) is the initial liquid phase contaminant concentration [g/ml], M (z,0) is the initial mass of contaminant per unit mass of soil at location z [g/g of soil], q is the volumetric water content, and r_b is the bulk density of soil [g/ml].

The concentration in the liquid phase is obtained by solving the following equation which accounts for advection given C_w(z,0):

In equation (2.5) is the Darcy's velocity of infiltration water. The gas phase concentration of contaminant is determined by the following diffusion equation:

where D is the effective diffusion coefficient of contaminant in gas phase. The above equation is solved with appropriate initial and boundary conditions. After the evaluation of C_w(z,t) and C_a(z,t), the equilibrium concentration of the contaminant in the three phases is determined as follows.

First the total mass in the soil M_T is calculated:

where f is the soil porosity. Then the individual phase concentrations are evaluated using the following relations.

The partial differential equations (2.5) and (2.6) are solved using an implicit, finite-difference numerical scheme. The distribution coefficient K_d [ml/g] and gas phase diffusion coefficient D [m²/day] are calculated using the following empirical relations:

where K_oc is the organic carbon coefficient of the contaminant and f_oc is the fraction organic carbon content of the soil, and D_air is the free-air diffusion coefficient of the contaminant.

Mathematical Discussion of VLEACHSM 1.0a

Lee (1996) added a term for the liquid-phase dispersion coefficient in the vadose zone, which is regarded as a linear function of the pore water velocity as:

where a is the longitudinal dispersivity [ft] of the vadose zone. Based on data from Gelhar et al. (1985) on the field measurements of dispersion in the unsaturated zone, Ünlü et al. (1992) used the following relationship in their model to describe the scale-dependent dispersivity:

where L_u is vertical distance from the source to the observation point ( i.e., the water table).

It should be recognized that dispersion remains a significant uncertainty in the modeling effort. By adding equation (2.8), (2.9), and (2.10) with (2.12), the one dimensional governing transport equation in VLEACHSM 1.0a is same as (Lee, 1996; USEPA, 1994):

where C_w is concentration of a contaminant in liquid (water) phase [mg/L]. C_a is concentration of a contaminant in vapor (air) phase [mg/L]. C_s is concentration of a contaminant in solid phase [mg/Kg]. q_w is volumetric water content (volume of water / total volume) [ft³/ft³]. q_a is air-filled porosity (volume of air / total volume) [ft³/ft³]. Note that the total porosity n equals the sum of water filled porosity and air filled porosity. q_w is water flow velocity (recharge rate) [ft/yr]. The air flow velocity (q_a) is assumed to be zero in this study. D_w is dispersion coefficient for the liquid phase contaminant in the pore water [ft²/yr]. D_a is gaseous phase diffusion coefficient in the pore air [ft²/yr]. m_w is first order decay rate of a contaminant in water phase [1/yr]. m_a is first order decay rate of a contaminant in gaseous phase [1/yr]. m_s is first order decay rate of a contaminant in solid phase [1/yr]. In this model, for simplicity, it is assumed that m_w = m_a = m_s = m. r_b is bulk density of the soil [gr/cm³]. z is vertical coordinate with positive being downward. t is time [yr].

2.2   Saturated Zone

2.2.1   Saturated Zone Mixing

The water table changes dynamically with time. As a result, the boundary between the vadose zone and saturated zone is dynamic. The problems of rising and falling water table require expanding or contracting the spatial domains of the models to accommodate the moving boundary. It is not a problem for the one-dimensional vadose zone model; however, for the saturated zone model it would require addition or deduction of nodes at the upper boundary. This would require an evaluation and switching of the set of nodes receiving fluxes from the vadose zone which appears undesirable. Dean and Carsel (1987) overlapped the spatial domains of two models and interpolate values for flux and transport of solute based upon the portion of water table. They said that it had the additional feature of eliminating the effects of the bottom boundary conditions prescribed for the vadose zone model on the simulation of solute transport just above the water table.

The saturated zone mixing process is simulated by employing a mass-balance principle (Lee, 1996). For the purpose of mixing calculation in the saturated zone, the conceptual vadose zone soil columns are assumed to be arranged in one of two ways. To represent a contaminant source in the soil, a group of vertical soil columns can be arranged at either transverse or parallel alignment to the groundwater flow direction (see Figure 2.3). For the transverse arrangement (Figure 2.3b), each soil column is treated independently both in the vadose and in the saturated zones. For the parallel case (Figure 2.3c), only the vadose zone migration in each column is treated independently. For the saturated zone mixing, for the parallel case, the up-gradient mixed concentration is used as influent concentration for estimating the next down-gradient column's mixing. In both cases, the bottom of each soil column is assumed to be parallel to the water table.

Figure 2.3 Schematic diagram, for parallel soil column arrangement case, illustrating the saturated zone mixing. The contaminant penetration depths underneath each soil column (Fig. 2.3b) were adjusted to produce a composite total depth (Fig. 2.3c) (Lee, 1996).

2.2.2   Groundwater Flow System

Flow nets are constructed using numerical solutions to the equations governing groundwater flow. Regardless of how a flow net is constructed, it provides not only a visualization of the groundwater flow paths, but also information which can be used to determine the rate of groundwater flow in a particular region.

Several computer codes (e.g., MODFLOW, FLOWPATH) are available which can incorporate two or three dimensional groundwater systems using detailed field-measured aquifer properties (McDonald and Harbaugh, 1988; Franz and Guiguer,1990; Hill, 1990; Carrera and Neuman, 1986a, 1986b, 1986c).

Toth (1962, 1963) found that the mathematical model is a realistic representation of the general configuration of the flow system where the topography is subdued and the water table slope is gentle. Toth (1963) also used a more general expression for the configuration of the water table in a region of gently rolling topography.

Governing Equation

The two-dimensional generalization of Darcy's law requires that the one-dimensional form be true for each of the x and y components of flow:

where h(x,y) is the head potential, q_x(x,y) and q_y(x,y) are the specific flow rates in the x and y direction, respectively, and K(x,y) is the hydraulic conductivity.

The two-dimensional continuity equation for steady-state conditions is

Laplace's equation combines Darcy's law and the continuity equation into a single second-order partial differential equation. Darcy's law, equations (2.17) and (2.18), is substituted component by component into equation (2.19) to give

If we assume hydraulic conductivity $K$ is constant, equation (2.20) becomes

Toth (1962) assumed that, in an undeveloped watershed, the fluctuations of the water table on an annual basis were small. That is, he used an average water table position and assumed that the system was at a steady state on an annual basis. The water table position at the beginning of the year was the same as the position at the end of the year; there was no net accumulation or loss of water from the system.

2.2.3   Advanced Groundwater Transport Models

Another main component of underground water system is incorporation of a saturated zone simulation model. The model should be able to handle both single (unconfined or confined) and leaky two-aquifer systems (Dean and Carsel, 1987). Transport processes of dissolved contaminants should include hydrodynamic dispersion, advection, linear equilibrium sorption, and first-order decay.

There is no doubt that groundwater movement is a three-dimensional problem. In three-dimensional cases, the governing partial differential equation is very complex and requires more conditions to solve it. It is often a common situation that the velocity of groundwater movement in vertical, z, direction is much smaller than that in horizontal, x and y, direction. Therefore, a most common assumption in solving groundwater problem is that there is no vertical flow in aquifer. Neglecting vertical movement of groundwater makes the solution of the governing equation much easier and keeps results accurate enough.

Therefore, an assumption that solute moves in two dimensions is made for the simulation in saturated zone. To solve the pollutant transportation in groundwater in two dimensions, three potential techniques could be considered. First, using conventional numerical technique to solve the governing solute transport partial differential equation; such as, Galerkin fate element method, Method of Characteristic (Koinkow and Bredehoeft, 1978), and Random Walk (Prickett et al., 1981). Second, using analytical element method (Strack, 1987; Strack et al, 1987), based upon superposition principle. It uses superposition to generate the analytical solution of a certain problem, expressed as a sum of basic solutions. The third way is to solve each subregion, which has homogeneous property, analytically then to integrate the obtained cell analytical solution to get a final solution. This technique is called cell analytical numerical method (CAN) (Elnaway and Valochi, 1989).

2.3   Model Linkage Problem

The connections of the vadose zone leaching, saturated zone mixing, and saturated zone flow sub-models in system is one of the most challenging problems. Both temporal and spatial linkages between the sub-models need to be considered, because unproper time step and grid space may result in divergent problem solution with the numerical technique.

2.3.1   Temporal Model Linkage

The time step for vadose zone model, for example, SAFTMOD (Huyakorn, et al., 1988) is dependent upon the properties of soils and the magnitude of the infiltration flux at the top of the soil column. The nonlinear Richard's equation might require time steps on the order of minutes to converge. In contrast, using finite element method to solve saturated zone model often requires much longer than one day time steps. In simulation, VLEACHSM 1.0a (Lee, 1996) the yearly time step is considered the basic time step, and correspondingly, using an integer multiple of yearly step as the time step of the finite different analysis in saturated zone. The daily infiltration output of PRZM (Carsel et al., 1984) is treated as a put input for vadose zone model at end of day. Similarly the relative continue output of vadose zone model is accumulated then put into the saturated zone model at the corresponding time step.

2.3.2   Spatial Model Linkage

The spatial linkages between models are more complex. The linking scheme should provide a realistic simulation of flow of water and transport of solutes at the interfaces and must ensure mass balance principle. The linkages between PRZM, vadose zone, and saturated zone models are different.

The VLEACHSM 1.0a uses simple "drainage rules" to move water through the soil profile. The major problem with interfacing of SAFTMOD and PRZM models is that Richards' equation is used for water flow in vadose zone and variably saturated medium. Because of this incompatibility, there may be times when root zone model produces too much out flow water for vadose zone model within one day. This is very likely to happen in agricultural soils, in which sub soil layers have typically lower permeability than those of the root zone are. In addition, the daily yield of PRZM might be unacceptable to the vadose zone model. The results of these would be water ponded at the interface which would belong neither to PRZM or vadose zone model. Dean and Carsel (1987) pointed out several options to handle this situation and prescribe the flux from PRZM into SAFTMOD so that the SAFTMOD accommodated all the water output by the root zone model each day. In this way, the problem of ponding at the two-zone interface was eliminated. But, on the other hand, it did force more water into the vadose zone than might actually occur in a real system. The consequence is that water and solute are forced to move at higher velocities in the upper portions of the vadose zone. If the vadose zone is deep, then this condition probably has little impact on the solution. If it is shallow, however, it could overestimate loading to groundwater, especially if chemical degradation rates are lower in the vadose zone than in the root zone.

The nodes linkage between vadose zone model and saturated zone model turn to the two dimensional (horizontal) model. Two models are linked through the use of a common node at the base of the aquifer, thus ensuring continuity of pressure head and mass balance of flux and solute.