3.1 Model Conceptualization, Assumptions, and Limitations .....3.1.1 Vadose Zone Leaching .....3.1.2 Saturated Zone Mixing .....3.1.3 Groundwater Flow System 3.2 Theoretical Background .....3.2.1 Vadose Zone Transport .....3.2.2 Saturated Zone Mixing .....3.2.3 Groundwater Flow System 3.3 Numerical Implementation .....3.3.1 Vadose Zone Leaching .....3.3.2 Saturated Zone Mixing .....3.3.3 Groundwater Flow System 3.4 Program Structure and Operation .....3.4.1 Input Parameters .....3.4.2 Input Data Format 3.5 Comparison between VG Model and VLEACHSM 1.0a .....3.5.1 Descriptions of Example Problem .....3.5.2 Results ## 3. MODEL DEVELOPMENT## 3.1 Model Conceptualization, Assumptions, and LimitationsAny mathematical model involves a simplification of real-world systems with the degree of simplification being inversely proportional to the sophistication of the model. Simpler models require more significant idealization of the real-world system. The Vadose-Groundwater Model (VG Model) is a relatively simple model intended primarily for screening contaminant level estimation within the solute transport phenomena in the subsurface. In this section, the assumptions and limitations are separately discussed for the vadose zone and saturated zone modules. ## 3.1.1 Vadose Zone LeachingThe 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 (see Figure 2.2). Equilibration among the phases occurs according to distribution coefficients defined by the user. Contaminant decay was incorporated in the model assuming a uniform decay rate for the three different phases of the user selected chemical. The first part (vadose zone leaching) of the VG model simulates vertical transport through vertically heterogeneous soil layers by advection and dispersion in the liquid phase and by gaseous diffusion in the vapor phase. These processes are conceptualized as occurring in a number of distinct, user-defined soil columns, which are vertically divided into a series of user-defined cells within each layer. The soil columns and layers may differ in soil properties, recharge rate, and depth to water. Within each soil column heterogeneous conditions are assumed due to differences in soil properties among different geologic layers and also in contaminant concentration (see Figure 3.1). During each time step intervals the migration of the contaminant through each cell in the soil column is calculated. Hence, the first part of the VG model can account for heterogeneity both laterally and vertically. The numbers of columns and layers employed in the simulation have a mathematical capability of 50 and 1000, respectively. Contaminant transport processes in the subsurface are simulated after
the first part of the
VG model has calculated contaminant equilibrium concentration in the liquid and
vapor phases based on the user provided initial solid-phase concentration [ For computational purposes each soil column is divided into a series of cells. When developing a model simulation, it is important to fully understand the implications of the first part of the VG model conceptualization. The following assumptions are made in the conceptualization of the first part of the VG model in the vadose zone. - Linear isotherms describe the partitioning of the pollutant among the liquid, vapor, and solid phases. Local or instantaneous equilibrium among these phases is assumed within each cell in every layer.
- Liquid phase dispersion and vapor phase diffusion are included. In the model the two terms are combined as one via the Henry's equilibrium constant.
- Unless
there is a special remedial action such as soil venting, vapor flow in the
subsurface is negligible. Although some
researchers (Falta
*et al.*, 1988) have found that density driven advective flow in the soil is also an important mechanism for gas movement, the vapor-phase advection is assumed as negligibly small compared to the infiltrating water velocity and omitted in this screening level model. - Organic contaminants, especially hydrocarbons, generally undergo some degree of degradation in the vadose zone. The module uses two decay constants to describe the first order decay; one for the contaminant decay in the soil, and another for the contaminant source reduction. The finite initial source concentration decreases due to either degradation or flushing by fresh water infiltration.
- There is no capillary fringe. In most permeable soils this is very small in proportion to the column length and can be neglected.
## 3.1.2 Saturated Zone MixingThe leachate from the vadose zone will reach groundwater in the saturated zone and will be mixed. In the second part (saturated zone mixing) of the VG model, this mixing process is simulated by employing a mass-balance principle. 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 Figures 3.2 and 3.3). For the transverse arrangement (Figure 3.2), each soil column is treated independently both in the vadose and in the saturated zones. For the parallel case (Figure 3.3), 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 concentration. In both cases, the bottom of each soil column is assumed to be parallel to the water table. For the screening level mixing calculation, the following assumptions are made in the development of the second part of the VG model. - The leachate is assumed to be completely mixed in the aquifer directly under the soil column.
- The aquifer is homogeneous with steady state and uniform groundwater flow.
- Groundwater mounding due to potential high infiltration is neglected.
- The decay of the solute in the saturated zone is conservatively assumed to be small and negligible.
- The leachate penetration in the aquifer is estimated using the relationship presented in EPACML (EPA, 1990a). It is introduced as a function of the vertical infiltration velocity, horizontal length of the source in the soil column, vertical transverse dispersivity in the aquifer, horizontal groundwater velocity, and the aquifer thickness. When the estimated penetration is greater than the thickness of the aquifer, the value is set as equal to the thickness of the aquifer.
- Contaminant spreading in the horizontal transverse
direction (
*i.e.*, perpendicular to the groundwater flow) is neglected. This assumption may cause a somewhat higher mixed concentration than would occur naturally.
## 3.1.3 Groundwater Flow SystemGroundwater flow models involving Laplace's equation and suitable boundary conditions are described by Toth (1962). He was able to draw conclusions about the configuration of regional groundwater systems by using a mathematical model. Figure 3.4 represents a cross section through a small watershed bounded on one side by a topographic high, which marks a regional groundwater divide, and on the other side by a major stream, which is a groundwater discharge area and marks another regional groundwater divide. The aquifer is assumed to consist of heterogeneous, isotropic, porous material underlain by impermeable rock. The
left and right groundwater divides can be represented mathematically as
impermeable, no-flow boundaries. Although no physical barrier exists, a
groundwater divide has the same effect as an impermeable barrier because no
groundwater crosses it. Groundwater to the right of the valley bottom
discharges at point A, and groundwater on either side of the topographic high
flows away from point In an undeveloped watershed, the fluctuations of the water table position can assume that the system is at a steady state on an annual basis, usually the U.S. Geological Survey water year (1 October to 30 September). The water table position at the beginning of the year is the same as the position at the end of the year; there was no net accumulation or loss of water from the system, viz, storage change.
## 3.2 Theoretical BackgroundIn this section, a brief theoretical development of the Vadose-Groundwater model (VG model) is presented for vadose zone transport, saturated zone mixing, and groundwater flow net.
## 3.2.1 Vadose Zone TransportIf we consider three equilibrium phases of pollutants in an unsaturated soil column, its one-dimensional heterogeneous governing transport equation is expressed as follows: 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_{w} is water flow velocity (recharge rate) [ft/yr].
The air flow velocity (q_{a}) is assumed to
be zero in this study. 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 is first
order decay rate of a contaminant in gaseous phase [1/_{a}yr]. m
is first
order decay rate of a contaminant in solid phase [1/_{s}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]. Instantaneous equilibrium (partitioning) of the contaminant among the phases is assumed to be realized according to the following linear relationships: Liquid-solid phase equilibrium is Liquid-gas phase equilibrium is where ml/g] is the distribution coefficient
between the solid phase and liquid
phase, and H (H = K /(_{H}RT), [dimensionless]) is the partition coefficient
between the air phase and water phase.
Using the empirical relationship, K can be expressed as _{d}K
= _{d}K·_{oc}f, where _{oc}K [_{oc}ml/g] is the organic carbon-water
partition coefficient and f [_{oc}g/g] is the fraction organic carbon of
the soil. The
dimensionless form of the Henry's partition coefficient, where atm-m^{3}/mol] is the dimensional form of Henry's Law constant,
R is
the universal gas constant (R = 8.2 · 10^{-5} [atm-m^{3}/mol·K]), and
T is the absolute temperature [°K =
273.16 + °C]. The
dispersion coefficient, The equations (2.14) and (2.15) can be rewritten for the English unit as: It should be recognized that dispersion remains a significant uncertainty in the modeling effort. Considerable debate exists regarding the appropriate practical treatment of field-scale dispersion. In general, field tests are the only methods by which reasonable estimates for dispersion can be obtained. The gas
phase diffusion coefficient ( where n
is the total porosity of the soil. By
substituting equation (3.2) and (3.3) into (3.1), and neglecting the air flow velocity
( If we set as By substituting equations (3.9) and (3.10) into (3.8) The equation (3.11) can be rewritten as: ## Initial and Boundary ConditionsTo solve equation (3.12), initial and boundary conditions are required as follows. where z,0) is the initial solid-phase concentration specified by the user.
When the distribution coefficient (K = _{d}K·
_{oc}f) is zero, liquid-phase concentration must be entered as
an initial concentration to avoid the program run-time error (division by
zero)._{oc}The most common type of boundary condition to be applied at the top of the soil column is either the first type (Dirichlet's) or the third type (Cauchy's) as shown below. or where At the bottom of the soil column, the second type boundary condition (Neuman's) is applied. In
applying this boundary condition, equation (3.19) is actually implemented at a
finite column length ( z,t), the total contaminant mass per unit volume of the
soil is calculated as:
## 3.2.2 Saturated Zone MixingAfter
estimating the liquid phase solute concentration ( et al., 1980):where q is the
Darcy velocity in the aquifer, _{aq}A is the cross-sectional aquifer area
perpendicular to the groundwater flow direction, and _{aq}A is the
cross-sectional
area perpendicular to the vertical infiltration in the soil column.
The aquifer area (_{soil}A) is determined by
multiplying the horizontal width of the soil column with the vertical solute
penetration depth._{aq}As discussed in Section 3.1.2, Saturated Zone Mixing, the procedure for the mixing calculation is different depending on the type of soil column arrangement. In the case of the transverse (right angle) arrangement, the mixing calculation is straight forward: simply apply equation (3.21) at each mixing element underneath the soil columns. For the parallel arrangement case, however, the mixed concentration at the upgradient cell is considered as an influx concentration to the next cell (see Figure 2.3). The mixing concentration at the next cell is estimated by re-applying equation (3.21) using the two inflow concentrations: one from the effluent from the upgradient cell and the other from the leaching from the soil column. ## Solute Penetration DepthThe solute penetration depth is the mixing thickness of the contaminant in the aquifer beneath the soil column. An estimation of the plume thickness in the aquifer can be made using the relationship below (adapted from USEPA, 1990a): where, ft],
a is the transverse
(vertical) dispersivity [_{v}ft] of the aquifer, L is the horizontal length
dimension of the waste [ft], and B is the aquifer thickness [ft].
In equation (3.22) the first term represents
the thickness of the plume due to vertical dispersion and the second term
represents that due to displacement from infiltration water.
When implementing this relationship, it is
necessary to specify that in the event the computed value of H
is greater
than B, the penetration thickness, _{d}H is set equal to _{d} B
(USEPA, 1990a).## Assumptions and Limitations Associated with the Penetration Depth EstimationThe
method introduced in equation (3.22) is based on two major assumptions (refer
to Appendix D for its derivation) and, therefore, has subsequent
limitations. First, in estimating the
dispersive penetration depth, there are two dispersion components in the mixing
zone: one from the transverse ( Secondly,
the second term in equation (3.22), the advective penetration depth due to
infiltration, is based on the assumption that the vertical (infiltration)
velocity varies linearly with aquifer depth, with a maximum value at the top of
the water table ( Due to
the assumptions and limitations discussed above, the user must be careful when
using the equation (3.22). Whenever
field data are available, the measured solute penetration depth should be used
instead of the equation. Otherwise, the
input parameter values must be checked to see whether the assumptions are
valid. For example, compare the
horizontal and vertical groundwater velocities to ensure the horizontal
velocity is much larger than the vertical one ( ## Adjustment of the Estimated Penetration DepthNo
adjustment of the penetration depth is required when either the depth is
measured in the field or the soil columns are arranged at right angle to the
horizontal groundwater flow direction.
For the parallel arrangement, however, when the penetration depth needs
to be estimated by the equation (3.22), the estimated depth is adjusted as follows.
As shown in the example (Fig. 2.3), the
total source dimension ( where
the subscript i represents column 1 through 4. Since the relationship between i.e., H
_{d }_{}
H + _{d1}H + _{d2}H
+ _{d3}H). Note that the overall
penetration depth (_{d4}H) is based on the total contaminant length dimension
(_{d}L) and the arithmetic average values of a,
_{v}B, q and _{aq}q for all columns. Therefore, the
individual penetration depth is adjusted as: _{w}where
Consequently, the relationship between the adjusted and individual penetration depth is expressed as:
## 3.2.3 Groundwater Flow System## Governing Equation and VisualizationThe
two-dimensional generalization of Darcy's law requires that the one-dimensional
form be true for each of the 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. Therefore, with this idealization, the two-dimensional Laplace equation (2.21) is the required governing equation. The mathematical model, consisting of the governing equation together with the four boundary conditions, is summarized in Figure 3.5. Because q satisfy equation (2.19) and thus (2.21),
they can be expressed
in the following forms:_{y}where For
visualization of the flow patterns, we have to define the flow rate vectors in
x- and y-directions. If h
function is given flow rate vector ( q) can be calculated from
equations (2.17) and (2.18). Then, _{y}V(x,y) can be calculated from equations (3.26) and (3.27).## Boundary ConditionsThe
coordinate system of Figure 3.4 is defined in
Figure 3.5. An equation is
required for each boundary. The upper boundary is located at The
other three boundary conditions are for no-flow boundaries. Darcy's law relates
flow to gradient of head. Along a vertical, no-flow boundary, dh/dx = 0, and along a horizontal, no-flow boundary, q
= 0 implies _{y}dh/dy = 0. Specification of flow across these
three boundaries makes them Neumann (second) boundary conditions. The full set
of boundary conditions is written as:The two
dimensional Laplace equation, Last modified: Oct 15, 1999 VG Model / Samuel Lee / VADOSE.NET |