2.1 Vadose Zone .....2.1.1 Structure of Vadose Zone .....2.1.2 Vadose Zone Transport .....2.1.3 Numerical Leaching Models 2.2 Saturated Zone .....2.2.1 Saturated Zone Mixing .....2.2.2 Groundwater Flow System .....2.2.3 Advanced Groundwater Transport Models 2.3 Model Linkage Problem .....2.3.1 Temporal Model Linkage .....2.3.2 Spatial Model Linkage ## 2. LITERATURE REVIEWA 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 ## 2.1 Vadose ZoneOrganic 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 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 ZoneThe 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,
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). ## 2.1.2 Vadose Zone Transport## LeachingLeaching 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 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
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). ## SorptionSorption 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 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, is the
soil organic carbon content (Wagenet and Rao, 1985
).
The large variability in _{f}K and
_{d}K due to adsorption by different soils can be greatly reduced by
normalizing adsorption based on the organic carbon content of the soil, _{f }K
(Jury, 1983)._{oc}## Spatial and Temporal VariabilitySpatial 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 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 ModelsVLEACH (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 VLEACHThe VLEACH treats spatial and temporal variation of contaminant concentration in
solid, liquid, and gas phases. These variables are denoted as 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 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 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 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 where The partial differential equations (2.5) and (2.6) are solved using an implicit,
finite-difference numerical scheme. The distribution coefficient ml/g]
and gas phase diffusion coefficient D [m^{2}/day] are calculated using the
following empirical relations:where f is
the fraction organic carbon content of the soil, and _{oc}D is the free-air
diffusion coefficient of the contaminant. _{air}## Mathematical Discussion of VLEACHSM 1.0aLee (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 where 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 mg/L]. C
is concentration of a contaminant in vapor (air) phase [_{a}mg/L]. C is
concentration
of a contaminant in solid phase [_{s}mg/Kg].
q is volumetric
water content (volume of water / total volume) [_{w}ft^{3}/ft^{3}].
q is
air-filled porosity (volume of air / total volume) [_{a}ft^{3}/ft^{3}]. Note that the
total porosity n equals the sum of water filled porosity and air filled
porosity. q is water flow velocity (recharge rate) [_{w}ft/yr]. The air flow velocity (q) is assumed to
be zero in this study. _{a}D is dispersion coefficient for the liquid phase
contaminant in the pore water [_{w}ft^{2}/yr].
D is gaseous phase diffusion coefficient in the pore air
[_{a}ft^{2}/yr]. m
is first order decay rate of a contaminant in water
phase [1/_{w}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 is bulk density of the
soil [_{b}gr/cm^{3}]. z is vertical coordinate with positive being downward.
t is time [yr].
## 2.2 Saturated Zone## 2.2.1 Saturated Zone MixingThe 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.
## 2.2.2 Groundwater Flow SystemSeveral computer codes ( 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 EquationThe two-dimensional generalization of Darcy's law requires that the one-dimensional
form be true for each of the where 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 ModelsAnother 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, 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 ## 2.3 Model Linkage ProblemThe 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 LinkageThe time step for vadose zone model, for example, SAFTMOD (Huyakorn, ## 2.3.2 Spatial Model LinkageThe 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. Last modified: Oct 15, 1999 VG Model/ Samuel Lee / VADOSE.NET |