CALIBRATED FLOW AND TRANSPORT MODEL FOR TCE IN GROUNDWATER AT BUILDING 460, FORT CHAFFEE, ARKANSAS

Prepared by: Darrel E. Dunn, Fort Lauderdale, Florida

Prepared for: Fort Chaffee Base Transition Team

August 10, 2001

[This is an html version of a report that was submitted in paper copy.  Transcribed by Darrel Dunn, Ph.D., Hydrogeologist, Colorado Springs, Colorado]

PURPOSE

This report describes the development of a groundwater flow and solute transport model at Building 460, Fort Chaffee, Arkansas, and the application of the model to calculating future concentrations of TCE at Little Vache Grasse Creek downgradient from the building.  The TCE was spilled in a northeast-southwest ditch along a road that passes northwest of the building.  The spills occurred sometime after the building was completed in 1942.

MODEL CODE SELECTION

The computer codes selected for the modeling were Modflow, Modpath, and MT3D.  These are widely known and accepted models.  Modflow and Modpath were developed by the U.S. Geological Survey.  MT3D was originally developed at S. S. Papadopulos & Associates, and subsequently documented for the Robert S. Kerr Environmental Research Laboratory of the U.S. Environmental Protection Agency.  MT3D was the preferred solute transport model because the most recent versions include the capability of simulating dual-domain porosity systems, in which immobile domain porosity contains immobile fluid where transport is primarily by molecular diffusion.  It was initially thought that it might be necessary to use dual-domain porosity to simulate the system.  This was not the case, and dual-domain porosity was not used.

GROUNDWATER FLOW MODEL

The groundwater flow model was developed and calibrated to supply groundwater flow velocities for particle tracking and solute transport modeling.

Conceptual Model

The subsurface materials in the vicinity of Building 460 consist of fractured shale overlain by clay.  A weathered zone has developed at the top of the shale.  The water table occurs in all three units at various locations around the site, but is commonly in the shale or the weathered zone.  The shale dips eastward at about one foot per forty feet and subcrops at an angular unconformity beneath the clay.  The dip of the shale is indicated by a dip-slope located northwest of Building 460.

The fractures may be conceptualized as a continuum of fractures with variable attributes.  The attributes that affect the groundwater flow include spacing between fractures, openness of the fractures, orientation, spatial continuity, wall roughness, and fracture fillings.  The depositional environment of shale changes from time to time, causing sedimentary layers with differences in composition and texture that affect the properties of the ultimate lithified shale.  These properties, in turn, affect the fracture attributes.  Consequently, some difference in fracture attributes may be expected to correspond to sedimentary layers, and the spatial distribution of fracture attributes will correlate with layering.  Furthermore, changes in depositional environment occurred laterally at any particular time, so that fracture attributes that affect the hydrologic properties of the shale may be expected to vary laterally within individual layers.  The fracture attributes that affect the hydrologic properties of the shale may be too subtle to identify by ordinary inspection and testing of cores and outcrops, although some, such as fracture spacing and orientation may be observed and described to a limited degree.  Rough estimates of hydraulic conductivity may be generated through in situ and laboratory testing.

The fracture attributes in the weathered zone have been altered during weathering, and the hydrologic properties of this layer may differ significantly from the unweathered shale.  The degree of weathering may vary from place to place due to changes in the parent rock fracturing and lithology.

The flow in the fractured rock and the weathered material derived from it is treated as conforming to Darcy's Law.  This is conventional practice because the equation relating flow velocity between parallel surfaces to separation of the surfaces has the form of the equation relating flow velocity through porous media to the hydraulic conductivity of the media.

Recharge to the groundwater system is through rainfall and snowmelt, and evapotranspiration returns water to the atmosphere.  Evapotranspiration extracts water from the land surface to the base of the zone affected by roots.  The water table conforms roughly to the topography, and the groundwater flows downslope toward Little Vache Grasse Creek.

The flow system is very complex.  Some observations at monitoring sites indicate increasing head with depth and others indicate decreasing head with depth.  Specifically, monitoring well MP03 penetrated a permeable, soft, broken shale at the bottom of the hole, which is 61 feet below the ground surface.  This is a multi-port well and the port in the deep permeable material has a head that is about 5 feet higher than that at the port about seven feet above it.  At MW05, head in the deep well is about 2 feet higher than in the adjacent shallow well.  Conversely, at the MW04 nest of wells, heads are more than a foot higher in a shallower well than an adjacent deeper one with only about 5 vertical feet separating the monitored intervals.  Such vertical head gradients require relatively low vertical hydraulic conductivities in the vicinity of the wells to sustain the head difference, and very heterogeneous lateral hydraulic conductivities.  The higher heads must have relatively high horizontal hydraulic conductivity present in the upgradient direction with low conductivity in the downgradient direction so that the head in the well is influenced more by the higher heads upgradient than by the lower heads downgradient.  The lower heads must have the opposite arrangement of lateral heterogeneity.  In all cases the flow is still downgradient to discharge to Little Vache Grasse Creek, but the head gradients are different in the upgradient and downgradient directions for the wells.  These opposing vertical head gradients can exist with the same upgradient and downgradient boundary  condition heads, but considerable heterogeneity in the hydraulic conductivity of the system is implied.

Discretization

A plan view of the spatial discretization of the groundwater model is shown in Figure 1.  The cell size is a uniform ten foot square, which was chosen to give good solute transport modeling results.  The boundaries were chosen far enough from Building 460 to make boundary effects insignificant in the area of interest.  The grid was oriented with the rows parallel to the ditch on the northwest side of Building 460.  This orientation produced lateral grid boundaries that are approximately parallel to flow lines, and allowed the TCE sources in the ditch to fall in a single row of model cells to facilitate the solute transport modeling.

The grid contains ten layers.  The top layer (number 1) represents the surficial clay.  Layers 2, 3, and 4 represent the weathered zone, and layers 5 through 10 represent the unweathered shale.  The weathered zone was divided into three layers to allow detailed simulation of the vertical distribution of TCE.  The bedrock was divided into six layers to allow detailed simulation of vertical head gradients measured in clustered and multiport monitoring wells.

Time was not discretized in the flow model, which was steady-state.

Boundary Conditions

The lateral and bottom boundaries of the model were set as no-flow boundaries, because they are approximately parallel to flow lines.

The upgradient boundary at the top of the model is a constant head boundary, with heads based on monitoring wells MW-03 and FTCLF2-MW03.  Linear interpolation was used between the monitoring wells.  Heads that were in cells not located between the monitoring wells are based on a water table elevations contour map of data from January and February 2000.  Constant heads in each column of cells on the top boundary are uniform (no change with depth).  Cells located above the water table were not assigned constant heads.

Little Vache Grasse Creek elevations were used as constant head elevations in the cells that are at the location and elevation of the creek.  The cells below the constant head cells were set as active cells, but the next cells to the south were set as no-flow cells.  This treatment at the creek makes the creek a dividing line between flow from the north and from the south.  Although no water level data is available south of the creek, such creeks normally function as sinks with flow directed toward the creek from both sides, flow beneath the creek being vertical.  If the flow is vertical, no horizontal flow exists and a vertical no-flow boundary may be placed at the stream.  The creek elevation passes from one lithologic unit to another along its course, so the constant head cells representing the creek are not all in the same model layer.  They are in the layer whose top and bottom elevations bracket the elevation of the creek at any particular place.

The treatment of Little Vache Grasse Creek as the downgradient discharge locus for the local groundwater system is consistent with known conditions at the site.  The Little Vache Grasse Creek valley is roughly symmetrical with the creek being the low point between two topographic divides that are about equal in height, roughly 80 feet above the stream level.  Monitoring well water level data show that the water table north of the creek follows the topography as is usually the case, rising to about 60 feet above the stream level.  The water table may be expected to follow the topography on the south side of the creek, as well.  Groundwater drains to topographic lows.  The resulting hydraulic gradient is toward the creek from both sides, and flow in the downgradient directions converges at the creek and must discharge to the creek.

The system being simulated is shallow.  The lowest cell bottom elevation in the model is 336 feet above sea level, which is about 76 feet below the stream level at that point, and 134 feet below the highest part of the water table on the north side of the creek.  The width of the valley between topographic divides is about one mile.  The flow system in the valley is 134 feet deep and a mile wide.  Such a shallow system should normally be controlled by local hydrologic conditions.

All water levels measured in wells north of the creek were at elevations higher than the creek, so the gradient at all depths in the shallow system under investigation is toward the creek.  One apparent exception to this observation is a low water level measured in MW06 adjacent to the stream.  The average water level measured in this well was 411.28, whereas the 412 topographic contour crosses the Little Vache Grasse Creek downstream from the well.  This discrepancy was interpreted as due to inaccurate stream elevations.  Contour elevations seem to be too inaccurate to use as exact hydraulic heads.  Another discrepancy between a measured elevation and one interpolated from the contours is at MW-8D which is given a land surface elevation of 427.14 but is next to the 426 contour.  This elevation problem at MW06 was dealt with in the model by setting the stream elevation next to MW06 slightly lower than the lowest water level elevation measured in the well.

Bottom Elevations

Bottom elevations of layers 1, 4,5, and 10 are shown in Figures 2 through 5.  The bottom elevations of layer 1 (clay) and layer 4 (weathered zone) were contoured based on borehole logs and elevations were stepped in two-foot increments.  Elevations of the bottoms of layers 2 and 3 were interpolated between the bottom elevations of layers 1 and 4.

The bottoms of layers 5 through 10 were derived by setting the bottom of layer 10 at 366 at MP03 and making the layer thicknesses agree with the port spacing in this well, one port per layer.  The layers strike north-south and dip eastward at one foot per forty feet.  Each layer has a uniform thickness except layer 5.  The bottom of layer 5 dips eastward like the other shale layers, but the top is the base of layer 4,which has an irregular elevation.  Consequently layer 5 becomes thicker toward the east.  Where layers are truncated beneath layer 4 in the western part of the model, they are extended to the western boundary as thin layers to accommodate the grid requirement that all layers are continuous through the grid.  These extensions did not affect the performance of the model.  They could have been divided into hydraulic conductivity zones equivalent to superjacent or subjacent hydraulic conductivities if necessary.

Where conflicts were found in layer bottom elevations (bottoms higher than tops), the elevations were adjusted the create layers that are at last two feet thick.

Hydraulic Conductivity

Hydraulic conductivity zones in the model were derived through calibration while keeping the values consistent with slug test results.  The preferential direction of the zones in layers 2, 3, 4, and 5 was north-south, following the strike of the shale because this trend would follow shale layers intersecting the ground surface or subcropping beneath the clay.  During calibration the north-south trend of the zones was preserved.  No directional trend for hydraulic conductivity zones in layer 1 or in layers 6 through 10 was imposed, because these model layers all follow the inclination of stratigraphic layers and no hydraulic conductivity trend can be inferred from the geologic conditions.  The hydraulic conductivity distribution of the calibrated model is shown in Figures 6 through 16.

Recharge

A uniform recharge rate of 0.000226 ft/day was used.  This is the estimated rate of percolation through the surficial material under unit head gradient, which is the gradient that corresponds to a condition of steady free vertical flow downward through the soil.  The value 0.000226 ft/day is the falling head permeability reported for a soil sample collected at a site north of Building 460 (U.S. Geological Survey, 1996, Section 3.2.1).  This value corresponds to 1.14 inches per year.  If the annual precipitation is 40 inches the recharge represents 2.85 percent of the precipitation.  If evaporation returns 85 to 95 percent of the precipitation to the atmosphere (U.S. Army Base Transition Team, 1998), then 0.86 to 4.86 inches runs off.

Calibration

Calibration of the model was accomplished by adjusting hydraulic conductivity in successive trials.  Calibration targets were average head at the monitoring wells listed in Table 1.  The calibration goals were as follows:

• Keep hydraulic conductivities within an order of magnitude of slug test results at the locations of the tests,

• Let hydraulic conductivities deviate as little as possible from slug test values to the north and south of each test in layers 3 through 5, where hydraulic conductivity trends might follow bedrock strike,

• Let hydraulic conductivities deviate as little as possible from the nearest slug test values in layers 1, and 6 through 10, and

• Maintain low flow across the upgradient boundary of the model, because the size of the recharge area upgradient from the model is small.

The chief performance measure applied to the calibration was residual standard deviation divided by the head range.  This value should be less than 10 to 15 percent for a good calibration (Rumbaugh and Rumbaugh, 1999, p. 257)  This statistic was 7.3 percent for the present calibrated model.  Other statistics for the calibrated flow model are given in Table 2.  The discrepancy in the volumetric water budget for the entire model was -0.74 percent.  Figures 6 to 16 show the hydraulic head, head residuals, and hydraulic conductivities for the calibrated model.

SOLUTE TRANSPORT MODEL

The results of the flow model were used to develop a solute transport model for the TCE plume at Building 460.

Discretization

The entire flow model grid was used for the solute transport model.  It was not changed.  The horizontal cell size remained at 10 feet by 10 feet.  The simulation was transient with a 58 year stress period representing the time from the beginning of 1942 to the end of 1999.  The conjugate gradient solver package was used, and transport time steps were calculated by the MT3DMS code.

Sources

The source cells are all in model row 62, which follows the location of the ditch northwest of Building 460.  Sources are point sources, simulated using input wells discharging into the source-cells at very low rates, so that the input rates would not impact the flow system.  The rate for each input well was 0.001 cubic feet per day.  Large input concentrations were used to compensate for the very low input flow rate and provide mass input rates that would not cause numerical problems.  The model input concentrations may be interpreted as grams per liter so that the calculated solute concentrations are grams per liter.  A multiplier of 106 converts to micrograms per liter (ug/L).  All calculated concentrations presented in this report are in ug/L.

Initially sources were placed in the model at single cells located where reverse particle tracking from the monitoring wells intersected row 62 or where vectors from the monitoring wells indicated the sources should be.  The sources were placed in the layer indicated by the particle tracking if the layer was 1 through 4.  If the reverse particle tracks intersected row 62 in layers below layer 4, the source was placed in layer 4,because vertical concentration profiles measured in the field showed concentrations peaking in the weathered zone.  The conceptual model is that most TCE below layer 4 reached that depth by vertical dispersion and diffusion from the weathered zone.  This conceptual model was used because concentration profiles measured in the field peaked in the weathered zone, suggesting that the TCE liquid generally did not penetrate into the shale.  Later in the calibration, some source cells were added in deeper layers because certain observed concentrations in the deeper layers required deeper source cells for calibration.  Source cells and their injection rates are shown in Figures 16 through 25.

Effective Porosity

The average measured porosity for the shale is 0.07.  The measured values do not vary much.  Consequently, a porosity of 0.07 was used in layers 5 through 10.  The porosity for layers 2, 3, and 4, representing the weathered zone was set at 0.10, because a porosity slightly greater than the bedrock porosity would be consistent with the DNAPL generally not entering the bedrock from the weathered zone.  The ten percent value is also the value recommended by Droppo and others (1989, p.234) for effective porosity of clay in the absence of measurements.  A value of 0.05 was used for layer 1, expecting the effective porosity of this compact clay the be less than the value for the fractured shale.

Diffusion

The diffusion coefficient for TCE is 9.4 x 106  cm2 /sec.  A medium tortuosity of 0.25 was used to obtain an effective diffusion coefficient of  2.35 x 106  cm2 /sec (2.185 ft2 /day).

Sorption

Initial distribution coefficients were calculated from laboratory values of fraction of organic carbon (foc) measured for subsurface samples collected at the Building 460 site.  The average values of foc for shale, weathered sale, and clay were 2.40, 1.85, and 0.525 ml/g respectively.  Using an organic carbon partition coefficient (Koc) of 126 ml/g, the calculated distribution coefficients (Kd) were 302, 233, and 66.1 ml/g.  These coefficients were initially applied to the appropriate layer.  However, Kd was extensively modified during calibration of the solute transport model.  The values were maintained as close to the initial values as possible.  The final values of Kd in the calibrated model are shown in Figures 16 through 25.

Biodegration

TCE is subject to biodegradation.  Decay half-life varies from site to site.  The present model uses a half-life of 300 days, which is given by Dragun (1988) as a value based on field observation in a naturally-occurring soil-groundwater system.  A test of the sensitivity of the model to the biodegradation decay coefficient was performed early in the calibration when only one source cell was active.  Values for half-life of 10-days, 300-days, 5-years, and 10-years were tested on this plume.  All but 10-days yielded virtually the same contours at the end of 58 years.  The value of 10-days shrank the plume greatly.

Calibration

The degree of calibration of the solute transport model may be seen on Figures 16 through 25.  The target concentrations shown on these maps are average values of TCE concentration for each monitoring well.  The contours are logarithmic, with the outer contour representing 1 ug/L.  The fit of the calculated values for TCE to the target values is within the range of uncertainty for the target values.  The uncertainty associated with the target values is large, as may be seen from successive values measured in the same monitoring well.  Values measured at MW01 were 1870, 9100, 9070, and 2630 ug/L.  At MW02 they were 230, 314, and 170.  At MW03 they were 942, 135, and 46.7.  The range can be more than an order of magnitude; and, since there are few measurements at a given site, the confidence interval for the true population mean is large.  This confidence interval was not calculated.

An important area in the calibration is between the ditch and nearest reach of Little Vache Grasse Creek, because this is the shortest contaminant travel distance.  In this area it was not possible to closely match target values, which ranged from 6 to 542 ug/L in layer 4, because the concentration pattern is not consistent with the normal distribution of concentration in a plume, which is declining concentrations away from the source and away from the plume axis.  In this case the plume contours were fitted to the target value of 125 ug/L at MP04.  The value 125 ug/L is virtually the same as the average of the values that could not be closely matched: 6, 13, 19, 51, 125, and 542 ug/L (the average being 126 ug/L).

Future Projection

The calibrated solute transport model was used to estimate future concentration of TCE in water entering Little Vache Grasse Creek.  To do this, all source concentrations were set to zero, the assumption being that the free TCE has disappeared from all source points or will be removed.  In that case all future evolution of the  plume is the migration of TCE already in the groundwater system toward the creek while subject to slow biodegradation.

Th simulated future plume moved very slowly.  Figure 26 shows the concentrations 32 years into the future in the layer containing the greatest concentration (layer 4).  Comparison with present concentrations shown in Figure 19 shows that the plume moved very little in the 32 years.  The higher concentrations decreased, and the 1 ug/L contour moved slightly toward Little Vache Grasse Creek.

Eventually, the plume reached the creek in the future runs.  The maximum impact on the creek was calculated to be about 1044 years into the future.  The location of maximum impact is at the meander near MW06, where the concentration reached 12.4 ug/L and then began to decline (Figure 27).  The concentrations in layer 5 at this time are shown in Figure 28.  The discharge of TCE into this reach of the stream at 1044 years is calculated to be 0.058 grams per day.

Another simulated TCE plume reached Little Vache Grasse Creek downstream near MW05.  The uppermost part of this plume is in layer 6 (Figure 29).  The creek is in layer 5 at this location and the discharge of TCE to the creek is upward from the cells beneath the boundary cells representing the creek.  The maximum TCE concentration beneath the creek was calculated at 29 ug/L.  However, the vertical hydraulic conductivity at this location is very low.  Consequently the discharge of TCE to this downstream reach is less than at the meander near MW06, even though the TCE concentration is greater.  The TCE discharge in this reach of the stream was calculated at 0.046 grams per day.  This plume is based solely on calibrating to concentrations measured in layers 6 through 9 in one well, MP05.

These two plumes account for all of the TCE reaching Little Vache Grasse Creek at 1044 years.

REFERENCES

Droppo, J. G. and others (1989): Multimedia Environmental Pollutant Assessment System Application Guidance, Volume 2, Guidelines for Evaluating MEPAS Input Parameters, Battelle Memorial Institute.

Rumbaugh, James O. and Douglas B. Rumbaugh (1999): Guide to Using Groundwater Vistas, Environmental Simulations, Inc.

U.S. Army Base Transition Team (1998): RI Report; Waste Accumulation Points (FTCH-021), Fort Chaffee, Arkansas.

U.S. Geological Survey (1996): Preliminary Report: RCRA Facility Invistigation, SMU-09, Former DDT Storage Site, U.S. Army Fort Chaffee Facility.

Walton, William C. (1985): Practical Aspects of Groundwater Modeling, William C Walton Consultant in Water Resources, 2nd Edition.