US Army Engineer R&D Center, Waterways Experiment Station,

3909 Halls Ferry Road, Vicksburg, Mississippi 39180-6199, USA

Abstract: An integrated modelling technique, which includes numerical simulation, optimization, and uncertainty analysis, was developed for investigating an optimal decision of contaminant plume capture design in groundwater. A generalized linkage between a numerical model, such as FEMWATER (a three-dimensional finite-element code) from the GMS (Department of Defense Groundwater Modelling System), and an optimization solver based on the unit-response matrix method was first constructed. The next step was to insert the parameter and boundary uncertainties into the hydraulic gradient control schemes. The developed sensitivity coefficients were then combined with the response coefficient calculation and considered as additional design constraints in the optimization process. A real case study to examine the optimal capture zone design for the uncertainties of hydraulic conductivity, infiltration rate, and constant head boundary was demonstrated.

Keywords: integrated modelling technique, uncertainties

Numerical simulation in groundwater has been commonly used as a decision-making tool for predicting flow and contamination transport and designing pump-and-treat systems. However, the available numerical tools are usually insufficient to design aquifer restoration needs and to meet specified management goals because of the tremendous number of trials that have to be repeatedly simulated. The combined simulation-optimization technique definitely provides quick decisions of well locations and pumping rates, however, previous investigations (Gorelick et al., 1984; Atwood & Gorelick, 1985; Colarullo et al., 1984), have all assumed that aquifer parameters are known with complete certainty, which never exists in any field study. To save project resources and to obtain a high probability of achieving project goals, an optimal design should include the integration of simulation, optimization, and uncertainty components.

Methods for the design of capture and containment systems can be developed using linear simulation management models if the hydrodynamic dispersion is ignored. Since there are always insufficient data on the physical parameters associated with pollutant movement through aquifers, another major objective is to decrease the sensitivity of the system to unforeseen changes in errors of the physical parameters. Valocchi & Eheart (1987) described a technique to incorporate parameter uncertainty into coupled linear optimization- groundwater simulation models.

In this paper, an extension using a three-dimensional numerical model and developing a generalized linkage between the numerical model and optimization schemes via hydraulic gradient control are reported. A real case study using with the investigations of the uncertainties of hydraulic conductivity in pumping the well layer, infiltration rate, and constant head boundary that affect the optimal design strategies was addressed.　

The system consists of three components, namely, numerical simulation module, optimization module, and uncertainty module.

The developed linkage between numerical simulation and optimization was a generalized procedure, which can fit both simulation tools in GMS. FEMWATER is a three-dimensional finite-element model of density-dependent flow and transport through saturated and unsaturated porous media.

The procedures for constructing a simulation model include the development of a hydrogeologic conceptual model, the development of computational mesh, the determination of boundary and initial conditions, the estimation of model parameters and performing the model run, the conducting of model calibration/verification, and the demonstration of pre-processing and post-processing.

The optimization module includes the steps of formulating the management problem, generating the response matrix, inputting the objective functions and constraints, solving the optimization problem, and post-processing of optimal results. This module uses the hydraulic gradient control method (Atwood & Gorelick, 1985; Gorelick et al., 1993), which is based on a unit-response matrix, to perform the linkage between the FEMWATER model and an optimization solver. The main idea for using hydraulic gradient control is to select wells and pumping or recharge rates so that the contaminant plume does not migrate during cleanup. As long as the plume is stationary, contaminated water can be removed, treated, and disposed or discharged. In the absence of hydraulic controls, the plume will migrate and disperse within the existing and future groundwater flow fields.

Minimizing well rates could mean to extract the most water or to inject the least water. Objectives based on well rates may include all or select managed wells. Once the management problem is converted into a mathematical formulation via response theory, it can be written in a format suitable for input into a linear or nonlinear programming solver. The original code, the revised simplex method, was adopted from numerical recipes (Press et al. 1988) and was extensively modified to conduct the construction of this module.

We assume that a best estimate of the aquifer parameter, such as hydraulic conductivity, has been in the traditional model described as the management scheme to determine an optimal pumping strategy. However, in some specific regions or layers, the hydraulic conductivity values are uncertain. We used the hydraulic gradient control concept to incorporate the sensitivity criterion into the management model by constraining the maximum possible value of the sensitivity coefficients. The sensitivity coefficients are used to measure the impact of a small change in this parameter upon contaminant plume containment.

The traditional management seeks to minimize the total injection and extraction rate while constraining the hydraulic gradient to be directed inward at the gradient check pair locations. The groundwater flow equation, such as the FEMWATER code, provides the coupling between the decision variable (Q_j) and the hydraulic gradients through the use of a response coefficient.

If we assume some specific regions of the aquifer parameter or the boundary values are uncertain, the sensitivity coefficients S_ik are used to measure the impact of a small change in this parameter upon contaminant plume containment. The sensitivity coefficient is defined as

S_ik = _¶ P_i / ¶ C_k                                                                            (1)

where P_i = h _Ai - h _Bi                                                                             (2)

is the hydraulic gradient at check pair i, and C_k is the parameter in layer k or boundary value.

The absolute value of S_ik is a measure of the robustness of any particular management strategy. That is, the strategy with small values of ABS(S_ik) will still successfully contain the contaminant plume even if the parameter or boundary value differs slightly from its assumed value.

Valocchi and Eheart (1987) incorporated the sensitivity criterion into the management model by constraining the maximum possible value of ABS(S_ik). The management model is given by

nw

Objective Function: Min S (u_j + v_j)                                                                             (3)

j=1

nw

Constraints: S R_ij (u_j -v_j) < g_i i =1, 2, … ngp                                                                            (4)

j=1

ABS(S_ik) < S_max                                                                             (5)

0 < u_{j <} Q_max j =1, 2, … nm                                                                             (6)

0 < v_{j <} Q_max j =1, 2, … nm                                                                            (7)

where u_j is the extraction rate at well j, and v_j is the injection rate at well j; g_i is the target hydraulic gradient at location i; R_ij is the gradient response coefficient; ngp is the number of gradient check points; nw is the total number of wells; Q_max is an upper bound on the maximum possible pumping rate; S_max is the absolute allowable sensitivity value. It is noted that the S_ik consists of a regional flow component and a gradient response coefficient component.

The developed system with three modules was applied to a restoration site to develop an optimal pumping strategy for determining the contaminant plume.

This remediation site is approximately 30 acres in size and contains several locations where past spills, disposal practices, and operations have contaminated soil and groundwater. This site and its surrounding vicinity is hydrologically characterized as a valley-fill aquifer system. The main aquifer appears to be unconfined and is observed to have a water table approximately 30 m below the ground surface at the site.

In response to TCE (trichloroethylene) contamination found in onsite wells, the study site installed a ground treatment facility and is continuing detailed studies of subsurface contamination at the site and its surroundings. Currently, groundwater extraction comes from five pumping wells with four of the five being deep-well pumps, with a daily pumping rate of about 1 million gallons. The treatment system is based on physical separation through air stripping to isolate TCE, with the treated water then being discharged into the river.

The numerical model used in this study was based on a hydrogeologic conceptual model. The extent of the esker deposit and geometry of the basement surface are of critical importance to the hydrogeological model. The map module of GMS was used to create a new surface 2D mesh (Fig. 1a). The basement and sediment geometry with nine different materials classified are shown with the 3D finite element mesh overlay in Fig. 1(b). Under the hydrogeologic condition modelled, the mesh is vertically divided into 9 materials and 21 layers. The 3D mesh consists of 11 000 nodes and 20 139 elements.

(a)                                                     (b)

Fig. 1  (a) 2D and   (b) 3D computational meshes.

The western and northwestern model boundaries are represented by the river and are defined by constant fixed head elevations of 116 m. The eastern boundary is a constant head of 165 m. The southern and northeastern boundaries are defined as no-flow boundaries. The recharge and water-well boundary conditions were also assigned to the model. A value of 1.32×10^-3 m day^-1 was used as a precipitation influx. Well-pumping boundary conditions were assigned to five locations. A numerical procedure was developed to estimate initial conditions with iteration process due to field measurements were not sufficiently dense enough.

Hydraulic conductivity was assumed to be homogeneous within each model layer. Unsaturated zone curves were required for moisture content, relative conductivity, and water capacity as a function of pressure head. The highest hydraulic conductivity is in the esker layer (86.4 m day^-1 for horizontal and 50.0 m day^-1 for vertical).

Since they’re presently is no planning for new wells and no injection activities for any pumping wells, hydraulic gradient control is based on the pumping strategies for these five existing wells. With the information from the past 2 years of pumping history, the average demand for each well and the total demand can be estimated.

The plume is located upgradient of the existing well system. Due to the hydrogeologic conditions and management strategies at the site, the hydraulic gradient control for preventing the plume migration was designed by assigning 15 control pairs on the down-gradient side of the well system. Each pair occupies two computational nodes; one is inward and other one is outward. The objective is to minimize pumping rates. Five constraints were associated with the capacity, demand, usage, and individual pumping rates and migration control. The constraints are as follows:

(1) Restrict the total volume pumped to allow a low-cost, low-capacity treatment system.

(2) Restrict the total volume pumped to meet cooling-system demands.

(3) Restrict the least usage of pumping well 3.

(4) Restrict individual pumping rates to allow the use of small, inexpensive pumps.

(5) Restrict heads outside the contaminated zone to values greater than or equal to the heads inside (no spreading).

To determine the response matrix, this study only requires six numerical simulation runs that significantly reduce the computer memory and time required obtaining the optimal solution. The optimal pumping solutions for these five wells are illustrated in Table 1. The total pumpage for satisfying the objective function and constraints is 7484 m³ day^-1. This amount is about 86 % of full-operation capacity of the entire treatment system. This preliminary result could indicate the current usage must be increased to meet all specified constraints.

It is necessary to verify the optimal solutions by checking the optimized values with respect to the desired constraints. Three different variables were used to present this process: velocity field, hydraulic head, and gradient at control pairs. Using the optimal pumping solutions as the pumping rate input and then rerunning the numerical model does this. The flow field and total head around the pumping and the hydraulic control area show that the resulting depression cone formed around the wells. Contaminated water will be trapped by these wells and will not be discharged to the river. The clearwater in the upgradient could bypass the contaminated area and discharge into the river. These results show that each control pair under the optimal pumping schemes satisfies the hydraulic gradient control strategy.

The sensitivity coefficient is defined as the total gradient change when the parameter or boundary value has some uncertainty associated with it. Actually, this total change is based on the regional hydraulic gradient change plus the total drawdown change due to the pumping activities.

Three experiments are investigated, 10 % (8.6 m day^-1) increase of hydraulic conductivity in the esker layer, 0.01×10^–3 m³ day^-1 of infiltration rate, and 1 m of constant head boundary over the northeaster side of the model increases. These changes are considered as the uncertainty due to the measurement or estimate error. For each experiment, six computer runs are made-one for background run due to parameter or boundary change and five for estimating the response coefficient matrix due to parameter or boundary change. Fifteen estimated sensitivity coefficients for each control pair were used as the sensitivity constraints for solving the optimization problem.

The S_max for the tests of each experiment needs to be estimated. Usually, the average response coefficient (S_avg) and the average drawdown can be used as the initially estimated S_max. The Table 2 summarizes the optimal pumping with given S_max for each experiment. Due to the scale difference among these parameter or boundary changes, the optimal pumping scheme changes are difficult to compare. However, for each individual parameter or boundary value, these results are very useful to determine what is allowable error in order to obtain reliable optimal pumping values. For those more sensitive control pairs, it might be the indication for more field measurement to support the reliable study. That some S_max obtain unfeasible optimal solution is due to the limitation of capacity constraint in the system.

A very effective and computationally efficient scheme was developed and then applied at a restoration site to determine the optimal pumping alternatives for groundwater remediation. An optimization procedure using GMS and FEMWATER was developed that could be used to contain the contaminant plume while minimizing pumping rates. The sensitivity constraint is used to examine the uncertainty of parameter or boundary value.

Considering the capacity, demand, individual pumping limit, usage priority, and migration control of contamination, the results indicate that the optimal total pumping (7484 m³ day^-1) was less than the full operation capacity (8646 m³ day^-1) at this demonstration site. However, from the optimal results with uncertainty constraint, it indicates that if the maximum sensitivity of hydraulic conductivity is less than 1.5×10^-3 or the maximum sensitivity of constant boundary head is less than 1.0×10^-2, total optimal pumping might need more than the water with full-operation capacity.

This work was funded by the U.S. Army Environmental Center. Permission to publish this information was granted by the Chief of Engineers.

Atwood, D., & Gorelick, S. (1985) Hydraulic Gradient Control for Groundwater Contaminant Removal. J. of Hydraul. 76, 85-106.

Colarullo, S. J., Heidari, M. & Maddock, III, T. (1984) Identification of an Optimal Groundwater Management Strategy in a Contaminated Aquifer Water Resour. Bull. 20(5), 747-760.

Gorelick, S., Freeze, A., Donohue, D. & Keely, J. (1993). Groundwater Contamination: Optimal Capture and Containment Lewis Publisher, FL.

Gorelick, S. M., Voss, I., Gill, P. E., Murry, W., Saunders, M. A. & Wright, M. H. (1984) Aquifer Reclamation Design: The Use of Contaminant Transport Simulation Combined with Nonlinear Programming Water Resour. Res. 20(4), 415-427.

Press, W., Flannery, B., Teukolsky, S. & Vetterling, W. (1988). Numerical Recipes Cambridge University Press.

Valocchi, A. J. & Eheart, J. W. (1987) Incorporating Parameter Uncertainty into Groundwater Quality Management Models. In: Systems Analysis in Water Quality Management (ed. By M. B. Beck), 251-257.