3.3 Numerical Implementation

3.3.1 Vadose Zone Leaching

The second derivative by depth z becomes

for 2<i<(n-1), and k = 1, 2, ���. 

If we define the dimensionless constants M_i, M_i', N_i, N_i', and L_i as:

equation (3.34) can be rewritten as:

Initial and Boundary Conditions

or

The finite difference forms of the top boundary conditions equations 3.15 through 3.18) for the soil column are: First Type Top Boundary Condition,

Third Type Top Boundary Condition (see Appendix C for derivation),

where

and M, N, and L were defined in equation (3.35).

The second type bottom boundary condition is used in this model as follows:

Note that as indicated in the section 3.2.1, equation (3.41) is supposed to be applied at the distance of z = . This finite column length will cause somewhat higher concentration estimation at longer times near the end of the column. To reduce this finite boundary effect, dummy cells can be added at the bottom of the column either by changing input file or the source code. When the source code is modified, estimated concentration at the water table (n-dummy cells) should be used in calculation of saturated zone mixing.

The finite difference forms of simultaneous equations (3.32 through 3.41) are computer coded in C++ to solve for the value of C_i^k+1 by the Thomas algorithm (Press et al., 1992).

Numerical Stability

Often, the efficiency of a numerical technique is limited due to the instability, oscillation, and mass-balance problems. Several methods have been proposed to determine the stability criteria of finite difference calculation (e.g., Fourier expansion method, matrix method, and other, [Hirsch, 1989]). The Fourier expansion method, developed by von Neumann, relies on a Fourier decomposition of the numerical solution in space neglecting boundary conditions. It provides necessary conditions for stability of constant coefficient problems regardless of the type of boundary condition (Mitchel and Griffiths, 1980). The matrix method, however, takes eigenvectors of the space-discretization operator, including the boundary conditions, as a basis for the representation of the spatial behavior of the solution (Ames, 1997; Hirsch, 1989). 

Based on the von Neumann method, Crank-Nicolson scheme of finite difference equation (Equations 3.28 through 3.34) can be derived as:

According to the stability criteria expressed in equations (3.42) and (3.43), it is clear that the combined dispersion coefficient (D) must be greater than zero. In natural soil conditions, it is rare to have a value of D as zero or close to zero. Specifically, if there is downward infiltration in the vadose zone, hydrodynamic dispersion of contaminant is inevitable. In addition, the air diffusion coefficients of selected organic compounds (in Appendix B) are must greater than zero. 

Note that the numerical criteria presented in equations (3.42) and (3.43) are general guidelines and may not work for all situations. To ensure a stable result, the output should be checked thoroughly. If the results show any oscillation or negative concentration, a smaller time step and/or cell size than the values from the above equations should be tried. After having "stable" results, additional runs are recommended to check convergence of the simulation results. In this case, even a smaller time step and/or cell size can be used to check whether the solution converges or not.

3.3.2 Saturated Zone Mixing

The saturated zone mixing equation (3.21) is computer coded in the second part (saturated zone mixing) of the VG program. The calculation results versus time will be showed in a graphic outputs on the second part of the VG program. 

Based on the mass balance principle of equation (3.21) and Figure 3.6, the mixed concentration is estimated as:

Figure 3.6 Illustration of saturated zone mixing mechanism for parallel column arrangement.

is the same as C_mx(i) in "Mixed Concentration in Groundwater" output data of the second part of the VG program. Similarly, mixed concentrations at the subsequent elements can be estimated easily.

3.3.3 Groundwater Flow System

Finite Differences Method

Numerical solutions yield values for only a predetermined, finite number of points in the problem domain. By limiting our need to know the head to a reasonable number of points N, we can convert a partial differential equation into a set of N algebraic equations involving N unknown potentials.

Figure 3.7 Finite difference grid of regional aquifer area, where the spacing in the x direction is Dx, the spacing in the y direction is Dy, and the aquifer thickness is the unit length.

Calculation of the Head Values (Gauss-Seidel Iteration)

A central approximation to d²h/dx² (x₀,y₀) is obtained by approximating the first derivative at [(x₀+Dx)/2,y₀] and at [x_0,(y₀-Dy)/2], and then obtaining the second derivative by taking a difference between the first derivatives at those points. That is,

which simplifies to

Similarly,

According to Laplace's equation (2.21), the finite difference approximation for Laplace's equation at the point (i,j) simplifies to

If the finite difference equation (3.49) were solved for h_i,j, then

where Dx and Dy are the length of x and y direction of cell. That is, the value of h_i,j at any point is an weighted average value of head computed from its four nearest neighbors in the nodal array. Equation (3.50) is often called the five-point operator because the algebraic equations that approximate Laplace's equation are created one after another by moving the star of five points throughout the domain of the problem. Suppose we guess a set of trial answers and then successively improve our guesses until we get the right answer. In short, iterative methods consist of guessing and adjusting. 

In Gauss-Seidel iteration method, we work through the grid in an orderly way; that is, we start at i = 0 and j = 0 and sweep across from left to right and down line-by-line as if we were reading a page. In this way, we can always use two newly computed values in the iteration formula. Thus the Gauss-Seidel iteration formula is

The change between two successive Gauss-Seidel iterations is called the residual c. The residual c is then defined by

By replacing h_i,j^m with h_i,j^m+1 after each calculation, the Gauss-Seidel procedure liquidates or relaxes the residuals at every node and thus leads to a solution of each algebraic equation. 

Let us consider the left boundary. The finite difference approximation to dh/dx is

For the left boundary, the point referred to by indices (i+1, j) is inside the problem domain, but the point referred to by indices (i-1, j) is outside. Therefore, we expand the finite difference problem domain by one additional column to the left by putting in a column of so-called imaginary nodes. The boundary condition dh/dx = 0 translates in finite difference form to h_i+1,j = h_i-1,j. That is, the value of head along the imaginary column must reflect across the left boundary. Thus, the left boundary is a symmetry line. The right and bottom no-flow boundaries are also handled by creating imaginary nodes. The top boundary has specified heads (Table 5.2) and, therefore, no imaginary nodes are needed along this boundary.

Calculation of the Flow Rates and Flow Lines

which is used to calculate flow rates q_x and q_y at each grid point from head values h of neighboring grid points. 

The finite differences expression of equations (3.26) and (3.27) becomes

In order to obtain flow function, V values from q_x and q_y based on equations (3.56) and (3.57), we use Gauss-Seidel iteration method as we did in this section to calculate the head values. The difference from this section is the boundary condition. Because we assume there is no-flow across left, right, and bottom boundaries, the flow function V should have a common value everywhere on those boundaries. Because that common value is arbitrary, we put V = 0 on those boundaries. Therefore, we initially put V = 0 at all the grid points and performed Gauss-Seidel iterations using equations (3.56) and (3.57) until residual error becomes below threshold. In the iteration process, the V values at the left, right, and bottom boundaries were kept 0, and the V values at all other grid points were optimized.

3.4 Program Structure and Operation

Figure 3.8 Flow chart of the Vadose-Groundwater (VG) model for computer programming.

3.4.1 Input Parameters

Simulation Data

1. NSCOL: Number of soil columns conceptualized for the site. Each soil column will have a unique set of parameter data. 
2. DELT: The model time step in years. By executing the preprocessor, user will be advised for selection of time step size based on the numerical stability criteria. 
3. STINE: The total length of the simulation time in years. 
4. TIMPL: Duration of contaminant source release in the unit of years. 
5. PTIME: The time interval at which the groundwater impact and mass balance results are showed to the "Mass Balance and Groundwater Impact Data" file in the output file screen. The output time interval is in years. 
6. PRTIME: The time interval at which the vertical concentration profile results are showed in the output file screen. The profile time interval is in years. 
7. ANGL: Angle of the soil column arrangement with respect to the saturated zone groundwater flow direction. A character either "T" or "P" should be entered to represent the "transverse angle" or "parallel angle" arrangement, respectively. 
8. KOC: The organic carbon distribution coefficient describes the partitioning of the contaminant with organic carbon. The coefficient is in units of [ml/g]. Appendix B lists the values of log(K_oc) for selected contaminants. If data regarding the pollutant being modeled is not presented, refer to the standard reference manuals that are documented in Appendix C or consult the manufacturer of the compound. 
9. KH: Henry's partition coefficient (H). H is an empirical constant that describes the liquid-gas partition of the contaminant. H is a function of the solubility and partial vapor pressure of the contaminant at a given temperature. The first part of the VG model utilizes the dimensionless form of H given as The dimensionless form of H can be determined from the more common form having the units of atmospheres-cubic meters per mole [atm-m³/mol] using the equation where H is dimensionless and K_H is in units of [atm-m³/mol]. Data regarding Henry's Law Constant for more than 60 common volatile and semi-volatile organic compounds are provided in Appendix B. 
10. CMAX: Values defining the maximum water solubility of the contaminant must have units of milligrams per liter [mg/L]. Appendix B provides water solubility information for more than 60 compounds. 
11. DAIR: The free air diffusion coefficient describes the transfer of the contaminant due to Brownian motion in the air phase [cm²/sec]. The preprocessor converts this metric unit to [ft²/yr]. Two tables in Appendix B list the values of the unlisted chemicals, refer to the standard handbooks such as Perry's Chemical Engineer's Handbook (Perry, 1984).

Soil Column Data

1. NCELL: The cell number defines the number of cells within the soil column (NCELL = Vertical Column Length / DELZ+1). 
2. DELZ: Vertical height of the cells within the soil column in feet. A guideline is provided in the preprocessor for selection of the cell height (Dz < 2D/q_w). 
3. WIDTH: Horizontal width dimension of the soil column in feet. 
4. LENGTH: Horizontal length of the soil column in feet. The horizontal area of the soil column equals to WIDTH times LENGTH. 
5. TYPTBC: Type of the top boundary condition of the soil column. In this model, user has a choice between the two types of conditions: 1 for the first type, or 3 for the third type. Note: For the bottom boundary, no option is allowed for users. The second type of boundary condition is used by the program. 
6. VALTBC: This parameter describes the value of the top boundary condition. For the first type boundary condition, a value representing constant liquid phase concentration at the first cell should be provided (C_w [mg/L] = C_s [mg/Kg] x 10^-6 / K_d [L/mg]). For the third type boundary condition, either the concentration of the infiltrating water or leachate source concentration (C_w [mg/L]) at the top of the soil column is need to be specified depending upon the field condition. 
7. ALDISP: The longitudinal dispersivity, a_v, in equations (3.5), in the vadose zone soil column [ft]. 
8. PLT : Variable to denote the output file option for soil and groundwater impact for plotting purpose. "Y" (or "y") indicates that the soil contaminant profile and groundwater impact will be shown in the output file screen. The vadose zone leaching contaminant concentration file contains the estimated contaminant concentration sorbed to the soil versus depth for specified time period. The mass rate of contaminant concentration file contains the mass rate of contaminant loading to the ground water versus time. 
9. PLTIME: Plot time defines the time in years for which the soil contaminant profile data will be created for the plot file. 
10. DECY1: Decay rate of the contaminant in the vadose zone soil column in the unit of [1/yr]. This parameter describes the degradation or natural decay of the pollutant. The decay rate (m) can be calculated from the contaminant's half-life using m = log_e(2)/t_1/2 = 0.693/t_1/2. Selection of this rate is not an easy task. In most cases, decay rates from the ideally controlled laboratory condition (Howard, et al., 1991; Sims et al., 1991) are higher than the rates measured from the field. Therefore, unless there is reliable field data available, a conservative value of zero is recommended. Otherwise, a sensitivity analysis for this parameter is necessary. 
11. DECY2: Reduction rate of the contaminant in the source [l/yr]. The finite initial contaminant source at the top of the soil column may decrease its concentration due to flushing by precipitation, degradation, or other factors. If the flushing is a main mechanism in source reduction, a simple mass balance technique can be applied to determine the value of DECY2. A conservative value of zero is recommended. Otherwise, a sensitivity analysis is necessary. 
12. J1 and J2: Upper and Lower cell number, respectively (J1 < J2). 
13. RECHRG: The groundwater recharge rate describes the velocity of water movement through the vadose zone. The rate is given in ft/year. In the vadose zone the hydraulic conductivity of the soil is a function of the water content of the soil. It should be noted that this parameter is extremely difficult to estimate as in reality it will vary with respect to time and space. It is strongly suggested that a range of possible recharge values be utilized to evaluate the potential variability of the results due to the uncertainty associated with this parameter. 
14. RHOB: The bulk density describes the mass of dry soil relative to the bulk volume of soil in units of grams per cubic centimeters [g/cm³]. Ranges for bulk density with respect to different soil types are given in Appendix B. 
15. POR: The effective porosity describes the volume of interconnected void space within the soil column through which water and air can flow. The effective porosity is less than the total porosity. Many rocks, crystal lines in particular, have high total porosity, but most of which may be unconnected. Effective porosity is a volume ratio in a dimensionless unit. 
16. THETA: The volumetric water content (water filled porosity) of the soil in percent total volume. This parameter is assumed constant in time and space, however, this rarely occurs in nature. The volumetric water content can neither exceed the porosity of the soil nor be lower than the irreducible (residual) soil water content [dimensionless]. 
17. FOC: The fraction of organic content of the soil is the relative amount of organic carbon present in the soil. This parameter defines the amount of potential absorptive sites for the contaminant in the solid phase and has a dimensionless unit. The fraction organic content can be determined from laboratory analyses or is documented in some soil classification and survey documents at a local Soil Conservation Service Agency Office. Generic values for organic content for soils of different texture are listed in Appendix B. 
18. XCON: Initial contaminant concentration in the column within a single or set of cells. The initial contaminant concentration must be defined for all cells within the soil column. Default is the solid-phase concentration [mg/Kg] assuming a non-zero value of distribution coefficient (K_d) is used (i.e., both KOC and FOC should be greater than zero). When K_d equals zero (i.e., KOC = 0 or FOC = 0), the sorbed concentration is always zero. In this case, liquid-phase concentration [mg/L] must be entered to avoid the numerical overflow (divide by zero) in the program. 

Aquifer-specific Data

1. FLXAQF: The Darcy velocity [ft/yr] of ground water (q_aq = volume flux per area = ki, where, k is the horizontal hydraulic conductivity and i is the hydraulic gradient). The area for the volume flux is the soil column width (WIDTH) times the contaminant penetration depth (PDEPTH).
2. CINAQF: The background concentration of the contaminant in groundwater in [mg/L].
3. PDEPTH: The contaminant penetration depth [ft] in the mixing zone. This parameter can be either measured from the field or estimated by the preprocessor using equation (2.20). Even if it is being estimated by the preprocessor, an arbitrary value (e.g., 1.0) must be provided. At the end of execution of this program, the estimated value of PDEPTH will be written over the arbitrary value.
4. PCHOIC: A character variable ("M" or "E") to specify whether the PDEPTH is measured or estimated. "M" stands for measured, and "E" stands for estimated. When PCHOIC equals to "M", the values of the next two variables can be arbitrary. However, when "E" is used for PCHOIC, appropriate values of the following two aquifer parameters should be provided to estimate the PDEPTH.
5. ADEPTH: The thickness of the aquifer in [ft]. The maximum vertical contaminant penetration depth in the aquifer is limited by the aquifer thickness. When the estimated penetration depth becomes greater than the aquifer thickness, the aquifer thickness is used as a contaminant mixing depth. Note that, even if the PDEPTH is measured, a numeric value must be provided for the execution of the preprocessor program.
6. AVDISP: The vertical transverse dispersivity, a_v, of the aquifer in [ft]. This parameter value is used for estimating the contaminant penetration depth by the preprocessor. As a rule of practice, the a_v is considered as 10 to 30 % of the longitudinal dispersivity, a_L. The a_L is usually estimated as to be 10 % of the travel distance (i.e., the horizontal length of the soil cell; Gelhar et al., 1985). Input an arbitrary number when PCHOIC equals "M". 

Groundwater Flow Sub-model Data File

3.4.2 Input Data Format

Each parameter value in the input data file must be contained within a field, which consists of 10 columns. In Table 3.2, an overview of the input data format is provided by listing each variable name in its proper record and field. As exceptions to a normal column width limitation of 10 characters, titles in records 1 and 5 can be as long as 80 characters. A complete description of the input data format is provided in Table 3.3 together with the definitions of the parameters. 

Table 3.2 Input Parameters for The VG Model

3.5 Comparison between VG Model and VLEACHSM 1.0a

3.5.1 Descriptions of Example Problem

Examples of the input data file (VLSM1A.INP and VLSM1B.INP) are provided in Table 3.3 and Table 3.4. The data input equivalent to VLSM1.INP file is reproduced as VLSM1A.INP in the VG model format. In VLSM1A.INP file, there is only one kind of layer with 51 pixels with bulk density r_b = 1.6, total porosity n = 0.4, and water-filled porosity q = 0.3. In VLSM1B.INP (see Table 3.4) file, the layer is divided into four soil layers whose total porosity decreases (from 0.44 to 0.38) and water-filled porosity increases (from 0.26 to 0.32) along with the depth shown in Figures 3.9 and 3.10. Figure 3.9 shows Case A where two soil columns are arranged transversely to the groundwater flow direction. Figure 3.10 shows Case B where two soil columns are arranged parallel to the groundwater flow direction. The contaminant concentration profiles in the vadose zone are the same between Case A and Case B, where only the saturated mixing zones are different. The bulk density is adjusted according to the total porosity change (contribution from the water content change is disregarded). The Soil Column 1 has 1st type top boundary condition. The Soil Column 2 (Case B) has 3rd type top boundary condition. This assumed set of parameters are derived from a geologic situation where the total porosity of soil decreases along with the depth due to gravitational pressure while the soil becomes wetter (water-filled porosity increases) along with the depth because water sinks down to the lower layers. Recharge rate q is kept constant (1.0 [ft³/yr/ft²]) in order to keep water-filled porosity of each layer constant. Organic content f_oc is also kept constant (0.005 [g/g]). 

Figure 3.9 Profile of Transverse Column Example Problem in the VG Model (where Soil Column 1 is the 1st type of boundary condition. Soil Column 2 is the 3rd type boundary condition).

Figure 3.10 Profile of Parallel Column Example Problem in the VG Model (where Soil Column 1 is the 1st type of boundary condition. Soil Column 2 is the 3rd type boundary condition).

 

3.5.2 Results

In both Soil Column 1 and Soil Column 2, the effects of putting four different layers instead of one vertically homogenous layer are obvious as shown in see Figures 3.11, 3.12, 3.13, and 3.14. Because the total porosity values for the first (0 ~ -5 [ft]) and second (-5 ~ -15 [ft]) layers are larger than the original value 0.4, liquid-phase contaminant went deeper after the same period of time (71 [mg/L] at -15 [ft] in Soil Column 1 at 10 years compared with 28 [mg/L] for the homogenous case). Because the total porosity and water-filled porosity values are the same with the homogenous ones for the third layer (-15 ~ -33 [ft]), the inclination of contaminant profile mostly seems identical to the homogenous case. After 20 years, the contaminant penetrates deeper while keeping a similar profile, and after 30 years, contaminant completely reaches the bottom of the vadose zone, mixing into groundwater.

Figure 3.11 Simulation Results of Vertically Homogeneous Soil Column I, 1st type of Boundary Condition at 10, 20, 30, and 40 years of Liquid-Phase Contaminant Concentration (mg/L) vs. Depth from Ground Surface (ft).

Figure 3.12 Simulation Results of Vertically Homogeneous Soil Column II, 3rd type of Boundary Condition at 10, 20, 30, and 40 years of Liquid-Phase Contaminant Concentration (mg/L) vs. Depth from Ground Surface (ft).

Figure 3.13 Simulation Results of Vertically Heterogeneous Soil Column I, 1st type of Boundary Condition at 10, 20, 30, and 40 years of Liquid-Phase Contaminant Concentration (mg/L) vs. Depth from Ground Surface (ft).

Figure 3.14 Simulation Results of Vertically Heterogeneous Soil Column II, 3rd type of Boundary Condition at 10, 20, 30, and 40 years of Liquid-Phase Contaminant Concentration (mg/L) vs. Depth from Ground Surface (ft).