Abstract
Background
The purpose of this work was to demonstrate an approach to groundwater remedial design that is automated, costeffective, and broadly applicable to contaminated aquifers in different geologic settings. The approach integrates modeling and optimization for use as a decision support framework for the optimal design of groundwater remediation systems employing pump and treat and reinjection technologies. The technology resulting from the implementation of the methodology, which we call PhysicsBased Management Optimization (PBMO), integrates physicsbased groundwater flow and transport models, management science, and nonlinear optimization tools to provide stakeholders with practical, optimized well placement locations and flow rates for remediating contaminated groundwater at complex sites.
Results
The algorithm implementation, verification, and effectiveness testing was conducted using groundwater conditions at the Umatilla Chemical Depot in Umatilla, Oregon, as a case study. This site was the subject of a governmentsponsored remedial optimization study. Our methodology identified the optimal solution 40 times faster than other methods, did not fail to perform when the physicsbased models failed to converge, and did not require human intervention during the solution search, in contrast to the other methods. The integration of the PBMO and Lipschitz Global Optimization (LGO) methods with standalone physically based models provides an approach that is applicable to a wide range of hydrogeological flow and transport settings.
Conclusions
The global optimization based solutions obtained from this study were similar to those found by others, providing method verification. Automation of the optimal search strategy combined with the reliability to overcome inherent difficulties of nonconvergence when using physics models in optimization promotes its usefulness. The application of our methodology to the Umatilla case study site represents a rigorous testing of our optimization methodology for handling groundwater remediation problems.
Keywords:
Groundwater modeling; Contamination; Remediation; OptimizationBackground
The increasing scarcity and degradation of potable water resources is an issue of global concern. Sources of water quality contamination include releases of chemical and radionuclide contaminants from pointsource and nonpoint source origins. The costs of addressing water quality issues on a global scale are substantial. Regulatory agencies such as the U.S. Environmental Protection Agency (EPA) and counterpart agencies in other countries must consider economic and human health costs with budgetary constraints in their efforts to clean up contaminated groundwater and restore water resources to beneficial reuse. The EPA documents this need to perform groundwater remediation in an optimal manner in the “National Strategy to Expand Superfund Optimization Practices from Site Assessment to Site Completion” (USEPA, 2012).
The objective of this research is to develop and demonstrate an automated, costeffective, and broadly applicable approach to groundwater remedial design applicable to contaminated aquifers in different geologic settings. This paper presents the development, scope and application of the simulationoptimization approach for remediating contaminated groundwater including reliability and efficiency verification to find globally optimized solutions to groundwater remediation problems. The approach is demonstrated using a wellstudied, publically documented site example that was the subject of a governmentsponsored remedial optimization study.
During the past two decades, researchers and engineers have been seeking methods for optimizing the design of ground water treatment systems. Table 1 provides an overview of existing algorithms and methods along with their salient features. Recognizing the limitations of the previously used optimization approaches and the need for efficient and effective groundwater remediation optimization tools motivates this research. This new approach combines flexible and efficient optimization techniques with commonly used subsurface groundwater flow and transport simulation models.
Table 1. Key features and limitations of previous optimization tools
DoD/ESTCP Simulationoptimization demonstration project
This approach is tested using a study problem posed as part of the joint U.S. Department of Defense (DoD)/Environmental Security Technology Certification Program (ESTCP) Groundwater Remediation Optimization Study (Minsker et al. 2004); the groundwater contamination remediation design at the Umatilla Chemical Depot in Umatilla, Oregon. The site has a preexisting and operational remedyinplace (RIP) installed to remediate the Royal Demolition Explosive (RDX) and 2,4,6Trinitrotoluene (TNT) contaminated groundwater plumes. The RIP consists of a groundwater pump and treat remediation system.
The original research teams used three optimization approaches during the DoD/ESTCP study. The DoD/ESTCP report presents these approaches in their entirety in Appendix D, Volume II. The design approaches consisted of Subject Matter Expertise (SME); the Modular Groundwater Optimization (MGO) (Zheng and Wang 2002), and; the Simulation/Optimization Modeling Optimization Software System (SOMOS) (Systems Simulation/Optimization Laboratory SSOL 2002). The team from the University of Alabama applied the MGO approach, the team from Utah State University applied the SOMOS approach, and the group from GeoTrans applied an SMEbased subjective engineering approach. The publically available project web site (http://www.frtr.gov/estcp webcite) provides the DoD/ESTCP study reports, groundwater flow and transport models, and modeling files of the final solutions to this problem.
The Umatilla research problem formulation used for testing is ESTCP Formulation 1. The goal of the formulation aims to reduce the projected clean up times of RDX and TNT at minimal cost, subject to constraints on the total allowable pumping and injection, treatment capacity, and the number of new wells needed. The experimental design of the DoD/ESTCP study directed the investigators to consider the existing flow and transport models as “uptodate and acceptable for design purposes”. The groundwater flow model code is USGS MODFLOW96 (McDonald and Harbaugh 1988; Harbaugh and McDonald 1996). The model code MT3DMS4 (Zheng and Wang 1999) for multispecies contaminant transport. The objective function calculator provided by ESTCP evaluates the cost associated with a remedial design and its performance. U.S. Army Corps of Engineers USACE (1996) developed and provided the MODFLOW96 groundwater flow and MT3DMS4 solute transport models to the ESTCP study group.
The algorithm developed in this study  PhysicsBased Management Optimization (PBMO) – is an automated simulationoptimization based method. Examination of the optimal design solution from PBMO with those from the ESTCP study demonstrates its effectiveness. The PBMO solution acceptance metric is a cost equal to or lower than developed by the MGO team. MGO provided one of the lowest costs with the least amount of computational effort of the automated approaches. Using the same physically based models and objective function calculator in all the studies isolates the performance of the demonstrated optimization algorithms.
Umatilla chemical depot, Umatilla Oregon site background and description
Briefly, Umatilla is a 19,728acre military reservation established in 1941 as an ordnance depot for the storage, renovation, and demilitarizing of conventional munitions, and for the storage of chemical munitions. As of 1994, the Umatilla site only stored chemical munitions awaiting destruction. A washout plant operated at the site in the 1950s and 1960s. Discharges to unlined lagoons consisted of an estimated 85 million gallons washout water laden with RDX and TNT. The water table is present about 47 feet beneath the bottoms of the lagoons. The resulting soil and groundwater contamination caused the placement of Umatilla on the EPA National Priorities List (NPL) in 1984. Section 3.1.1 of the study report (Minsker et al. 2004) provides additional description of the Umatilla site and historical summary.
The Record of Decision (U.S. Army Corps of Engineers USACE 1994) – the legal document governing the remediation approach  specified a pump and treat system with reinjection of treated groundwater as the remedial alternative. The treatment system commenced operations in January 1997; this action is the RIP. The system comprised three active extraction wells [EW1, EW3 and EW4], three active infiltration/recharge basins [IF1, IF2 and IF3], and a granular activated carbon (GAC) treatment system with a capacity of 1,300 gallons per minute (gpm). Figure 1 provides the locations of the system components. The extraction well (EW2) installed 100 feet northwest of EW4 is not part of the RIP, but is available for inclusion as merited.
Figure 1. Initials conditions at end of 2002 extraction well and infiltration.
The initial RIP included an existing industrial lagoon designated IFL. However, suspicions that infiltrating water could spread the TNT plume rather than flushing the unsaturated zone of contamination resulted in discontinuing its use. By midJuly 1999, the RIP system had extracted, treated and recharged approximately 1.27 billion gallons of groundwater and removed an estimated 3,000 kilograms (kg) (6,614 pounds) of RDX and 400 kg (882 pounds) of TNT. The simulated RDX and TNT contaminant plumes, shown in Figure 2, represent the extent of contamination in the shallow aquifer at the end of 2002. The maximum RDX and TNT concentrations at that time were 28.2 and 86.7 micrograms per liter (μg/L), respectively.
Figure 2. Initials conditions at end of 2002 new well search regions candidate infliction basins and hydrogeologic yield zones defined for MGO.
Methods
Strategic planning for groundwater remediation, water resources planning and dewatering for mineral and resource mining requires tools that incorporate the complex constraints associated with the environment and support decision making with varying levels of uncertainty. This can include uncertainty in the conceptual site model, subsurface characterization (e.g. geologic material location and properties), chemical transport (e.g. reactions, natural attenuation, biodegradation), and funding levels (e.g. annual and total life cycle project) (ITRC, 2007).
Physicsbased models: MODFLOW, MT3D, MODFLOWSURFACT, and; MODHMS (HydroGeoLogic, Inc. HGL 2012) are highly effective at integrating remedial technology selection and application with competing remediation goals among regulators, site owners/custodians, and other stakeholders. The reasons for this are:
1.
 Comprehensive
2.
 Efficient and Effective
3.
 Flexible
Groundwater remediation problems amenable to physicsbased decision optimization include the development of costeffective and sustainable remedial designs; RIP evaluation and operation and maintenance (O&M) optimization; the development of optimized exit strategies to minimize lifecycle costs; the development of site completion/closure strategies; and water quality management issues for riparian and lacustrine settings, wetlands, and estuaries. This approach provides decision support to program and project managers regarding the best methods for remediating a site with contaminated groundwater.
Cost minimization is the overall optimization objective. In the case of groundwater remediation, well locations and their extraction or injection rates are the decision variables, and the simulated hydraulic heads and contaminant concentrations in groundwater are the state variables. Decision variables are the design elements of the problem studied. State variables represent the resultant simulated system. The objective function is a systematic accounting of costs which include the number, operation, and maintenance of the design elements (decision variables) and account for the changes in the modeled system (state variables) response from the candidate design. Computed over the life cycle of the simulation, the objective function represents the total remediation cost. The incorporation of constraints ensures the solution is practical and implementable. Constraints capture mechanical or physical process limitations (such as maximum pumping rates, and treatment train capacities) or interim and final values of state variables at discrete or continuous regions of the modeled domain. Assessment of candidate solutions is conducted by examining the cost value and the constraint compliance. A single model function evaluation is the cycle of the decision variable selection, process simulation, model objective function evaluation and constraint evaluation of a candidate design. Optimization is the automatic, guided process of the decision variable selection and application through repeated model function evaluations to arrive at the optimal objective function value that satisfies all constraints. The automatic adjustment of decisions variables terminates (ideally) at the global optimum (least cost).
PBMO’s approach provides decision support for groundwater remediation, water resources planning and dewatering for mineral and resource mining via flexible integration of physicsbased (groundwater flow and transport) models and optimization algorithms. Modular design promotes flexibility via independence of the physics model and optimization method. Linking the appropriate physicsbased simulator with the best optimization algorithm(s) enables the solution to a wide range of problems types. Figure 3 shows how the modular design links the global optimization algorithms with the appropriate physicsbased simulators to develop an optimal strategy. In this work, the first (top) half of the medallion represents the HGL_OPT optimization algorithms which include the Lipschitz Global Optimization (LGO^{©}) solver suite (Pintér 1996 2002 2009 2013). The second (bottom) half represents the physically based numerical groundwater flow and transport simulators USGS MODFLOW96 and MT3DMS4B (in this study); however any process simulator with a text interface is implementable.
Figure 3. PBMO process flow diagram.
Subsurface simulators represent the physical processes of flow and transport in a mathematical form. Problemspecific, physicsbased models represent the groundwater processes/conditions under investigation. From an optimization perspective, these simulation models can be “blackboxes”, that is models that receive input to produce an output without revealing the process. Since these groundwater flow and transport processes can occur in saturated or variably saturated conditions, in porous or fractured media, the numerical model used for the simulation will have a linearity category; the simulation can be linear, mildly nonlinear (assumedly or provably convex), or highly nonlinear (assumedly or provably nonconvex). The system under consideration can be either single phase (water) or multiphase (waterNAPLgas). The transport processes can be represented simply by using advective particle tracking, as a single nonreactive solute that undergoes dispersion and diffusion, or as a multicomponent system whose constituents interact with each other and the surrounding environment and are subject to hysteresis. The equations that govern these processes are Lipschitzcontinuous, whose statespace representation can range from elliptic to parabolic or nearly hyperbolic when advection dominated. When no or low diffusion and dispersion are present relative to advection, steep fronts or shocks occur in the solutions. A complication for efficient optimization approach design arises when that the decision variables used in a candidate model evaluation during optimization change the model linearity category. Typically, the simulation of groundwater flow with slight drawdown is linear to mildly nonlinear (because the saturated aquifer thickness does not change appreciably), whereas the simulation of significant groundwater drawdown or coupled groundwater flow and transport are most often highly nonlinear. For the model scenario considered in this study, the flow system will range from linear far from the remediation system and contaminant plumes to highly nonlinear in the contaminant plumes and near the remediation wells. The objective function is nonlinear.
The USGS numerical simulator, MODFLOW96, solves the following 3D governing partialdifferential equation for transient, saturated groundwater movement through porous media to yield a spatial distribution of the potentiometric head as a function of time:
Where K_{xx}, K_{yy}and K_{zz} are the values of hydraulic conductivity [L/T] along the x, y, and z coordinate axes, respectively, assumed to be parallel to the major axes of hydraulic conductivity; h is the potentiometric head [L]; q_{s} is a volumetric flux per unit volume of the aquifer and represents sources (injection wells) and/or sinks (extraction wells) of water [1/T]; Ss is the specific storage of the aquifer material [1/L]; and t is time (T).
Injection wells represent the infiltration/recharge basins. The initial ESTCP study used some terms interchangeably: injection, infiltration, and recharge. For clarity with the initial study language, we maintain that terminology in this paper to refer to the return of treated groundwater to the subsurface. The distribution of heads (one of the state variables), the source/sink data and, the candidate decision variables used or generated from the MODFLOW96 simulation provide input to MT3DMS to solve the partial differential equation describing the fate and transport of contaminant species k in 3D, transient groundwater flow systems. The governing equation is:
Where θ is the porosity of the aquifer material []; C^{k} is the dissolved concentration of species k [M/L^{3}]; x_{i,j} is the distance along the respective coordinate axis, x, y or z [L]; D_{ij} is the hydrodynamic dispersion coefficient tensor [L^{2}/T]; v_{i} is the seepage or linear pore water velocity [L/T]; is the concentration of the source or sink flux for species k [M/L^{3}], and; ∑ R_{n} is the chemical reaction term [M/(L^{3}T)].
Darcy’s Law couples the transport and flow equations.
McDonald and Harbaugh (1988) provide details on MODFLOW. Zheng and Wang (1999) provide details for MT3DMS.
To provide an optimization solver with the capability and efficiency to address this range of model types, PBMO is developed and implemented leveraging a suite of optimization algorithms that efficiently solve the various formulations expected to occur in the optimal decision support approaches presented here.
HGL_OPT is the master optimization driver. It comprises machine learning, objective function estimating methods, and SMEdeveloped heuristics specifically useful for quickly developing high quality solutions to complex problems such as contaminant flow and transport remedial design (Deschaine 1992 and Deschaine 2003; Deschaine and Pintér 2003; and Deschaine et al. 1998 2001 and Deschaine et al. 2011). These solutions can be used to focus (or initialize) the LGOspecific global and local optimization solvers.
The optimization algorithms included in the optimizer suite include linear programming (LP), sequential linear programming (SLP), sequential linear approximation (SLA), sequential quadratic programming (SQP), generalized reduced gradient (GRG), outer approximation (OA), Branch and Bound (BB), Globally Adaptive Random Search (GARS), and Multistart Random Search (MS). Specifically, LP solves linear models; OA, SLP, SLA, SQP, and GRG serve to deal with mildly nonlinear models; BB, GARS, and MS―in proper combination with SLA, SQP, and GRG―serve to handle highly nonlinear models.
Currently, the optimization system can handle model formulations with up to a few hundred binary (yes/no) decision variables, several thousand continuous variables, and a several thousand general constraints (in addition to variable bound constraints). The actual configuration can be adjusted to user demands. For example, it is possible to optimize models with 200 decision variables and 5,000 general constraints, or vice versa. Even larger model sizes can also be handled if the computer randomaccess memory (RAM) and the compiler used support this. The BB solver of LGO employs a mild assumption of Lipschitz continuity, which allows an efficient and a robust search for a global optimal value through a systematic partitioning and exploration of the entire feasible region. The GARS and MS solvers assume basic model function continuity. The optimization suite solves the following general model:
i. Minimize the [arbitrary linear or nonlinear] objective function f(x)
ii. Subject to bound constraints and [arbitrary linear or nonlinear] constraints:
where x is a decision vector, an element of the real nspace R^{n}; f(x) is a continuous objective function, f:R^{n}R, (R=R^{1}); D is a nonempty set of admissible decisions, a subset of R^{n}. As shown by (4), the set D is defined by l, u which is an explicit, finite nvector bound of x (a “box”) in R^{n}; and g(x) which is an mvector of additional continuous constraint functions, g:R^{n}R^{m}.
The assumption is that Lipschitz continuity holds as expressed by the inequality
where [L_{j}] is the Lipschitz constant. This equation means that the variability of the values computed by the objective function calculator is bounded with respect to the variability of the system input variables [x] by the Lipschitz constant. This condition is an expected property of all physically based models.
Effective development and specification of a problem for groundwater remediation projects requires SME understanding of the physical processes and their constraints. The options under consideration include selecting which of the physicsbased processes that need modeled (such as aquifer quality remediation via bioremediation or pumpandtreat) along with the constraints that the solution must satisfy. The formulation below is directly applicable and extensible to a wide range of groundwater dewatering and remediation problems. The basic elements of the mathematical formulation consist of an objective function and constraints.
subject to:
Here, the objective function f(x) is a simplified version of life cycle cost for illustrative purposes. It consists of the flow rate from extraction wells [q_{i}], a unit cost to process and treat the water [α_{i}]; and it assesses if an extraction well is preexisting or not by examining for ∀ i locations and assigning a (0,1) multiplier if a well installation required (δ_{i}=1) or is preexisting (δ_{i}=0). The actual formulations of the objective function for optimal design problems often involve quite extensive and detailed cost information and functions such as used in this study. The constraint requires the installed system pumps at least some minimum volume of water Q. The constraint sets upper limits on a flow at each of the i^{th} candidate well. The constraints ensure aquifer remediation to a certain acceptable residual level for all j locations. The constraints set minimum allowable water table elevations in the aquifer at all k locations. The t ≤ T_{max} constraint limits the allowable time for the activity to achieve the specified goals. These general constraints, like the cost function, can be modified and adapted as dictated by the needs of the project. For example, one can use this framework to incorporate constraints for differential land subsidence due to dewatering by using differencing constraints and adding constraints on the maximum slope of the water table or land surface. Similar approaches are effective for assessing depressurization and geotechnical stability.
The goal of constrained global optimization is to find at least one point x^{*} within the feasible region that satisfies f(x^{*}) ≤ f(x) for all x or to show that such a point does not exist. If no solution exists, then the decision makers realize it is an infeasible problem formulation. This enables the problem be reformulated. It is critical to determine if a feasible solution does not exist when solving a complex decision problem. This is a valuable capability of the physicsbased optimization methodology. Many other optimization methods, specifically including heuristic techniques, cannot determine that a problem is infeasible: as a result, users can waste precious time and money searching for solutions with no hope of finding a feasible solution.
Umatilla optimization problem formulation
The three optimization problem formulations considered in the DoD/ESTCP study were:
Formulation 1: Minimize the cost to remediate RDX and TNT in 20 years or less using the current treatment plant maximum operating flow rate of 1,300 (gpm) as an upper limit on the total groundwater extraction rate;
Formulation 2: Same as Formulation 1, except the maximum treatment flow rate increased to 1,950 (gpm);
Formulation 3: Minimize the aggregate remaining mass of RDX and TNT in 20 years using the current treatment plant flow rate.
The PBMO approach addresses all these types of optimization formulations. This exercise focuses on finding a solution to the problem as stated by Formulation 1. Appendix D, Volume II of the DoD/ESTCP report (Minsker et al. 2004) provides the mathematical formulation of the problem statement in detail.
The objective function for Formulation 1 is specified as follows:
where:
C_{CW}: Capital costs of new wells
N_{y} is the modeling year when cleanup is achieved [yr] as defined by [C_{RDX}] ≤ 2.1 μg/L and [C_{TNT}] ≤ 2.8 μg/L as measured by the nodal concentration value in the top layer.
N_{Wi} is the total number of new extraction wells (except well EW2) installed in year i. New wells may only be installed in years corresponding to the beginning of a 5year management period. Capital costs do not apply to preexisting extraction wells.
I_{EW2} is a flag indicator; 1 when well EW2 first comes into service, 0 otherwise.
75 is the cost of installing a new well [K$].
25 is the cost of putting existing well EW2 into service [K$].
d indicates the application of the discount function to yield Net Present Value (N_{PV}) defined as
c is the cost
r is the annual discount rate [1/yr]
y is the value of i in the summation
C_{CB}: Capital costs of new infiltration basins
N_{Bi} is the total number of new infiltration basins installed in year i. New recharge basins may only be installed in years corresponding to the beginning of a 5year management period. The infiltration flux is evenly distributed throughout the basin.
25 is the cost of installing a new recharge basin independent of its location [K$/yr].
F_{CL}: Fixed cost of labor
237 is the fixed annual O&M labor cost [K$/yr].
F_{CE}: Fixed cost of electricity (lighting, heating, and the like).
3.6 is the fixed annual electric cost [K$/yr].
V_{CE}: Variable electric costs of operating wells (extraction and injection)
is the total number of extraction/injection wells active in year i.
C_{Wij} is the electrical cost for well j in year i. Costs differ for wells depending on the extraction rates of well j in year i, Q_{ij}:
I_{Wij} is a flag indicator; 1 if well j is active in year i, 0 otherwise.
V_{CG}: Variable costs of changing GAC units
is the average influent concentration [μg/L] (RDX plus TNT) into the treatment plant from all of the extraction wells in the year i calculated as:
is the average influent concentration [μg/L] (RDX plus TNT) from well j in the year i
is the cost of mass removed [K$/kg] as a function of average influent concentration into the treatment plant in year i, calculated as:
m_{i} is the mass of contaminant removed [kg] during year i calculated as:
β is a conversion factor to produce the result in units of [kg/yr].
V_{CS}: Variable cost of sampling.
I_{A} is the initial plume area in layer 1 of the model based on the extent of RDX and TNT as measured in January 2003 where RDX and TNT exceeded their respective cleanup goals (2.1 and 2.8 μg/L, respectively) [m^{2}]
150 is the annual sampling cost (as of January 2001) and considers both labor and analysis costs [K$/yr]
A_{i} is the modeled plume area in layer 1 in year i. This is evaluated at the beginning of a 5 year management period [m^{2}]. A_{i} is defined as:
N_{col} is the number of grid cell columns in the x direction
N_{row} is the number of grid cell rows in the y direction
Δx_{j} is the width of the j^{th} grid cell column [m]
Δy_{k} is the width of the k^{th} grid cell row [m]
is the concentration of RDX in the grid cell with indices j and k
is the concentration of TNT in the grid cell with indices j and k
The Formulation 1 constraints are:
1) The modeling period consists of four 5year management periods beginning with January 2003 (i or year = 1).
2) Modifications to the system may only occur at the beginning of each management period.
3) Remediation in the top layer of the model must be achieved within 20 years (e.g., RDX ≤ 2.1 μg/L and TNT ≤ 2.8 μg/L everywhere in top model layer).
4) The total pumping rate, adjusted for the average amount of uptime, cannot exceed the treatment capacity of 1,300 (gpm) in any stress period. Evaluation of this constraint occurs at the beginning of each 5year management period. It is computed as:
Here
α is a coefficient that accounts for the amount of average uptime (α=0.9)
Q* is the total modeled flow rate during a 5year management period.
5) The hydrology dictates the upper sustainable flow limit on extraction wells. Extraction wells in Zone 1 may pump at a maximum rate of 400 (gpm), whereas extraction wells in Zone 2 may operate to a maximum of 1,000 (gpm). See Figure 2 for definitions of Zones 1 and 2:
where:
Zone(j,k) is a function of the j^{th} grid cell column and k^{th} grid cell row that returns 1 if model grid (j,k) corresponds to Zone 1, and returns 2 if (j,k) corresponds to Zone 2
is the modeled extraction rate at model grid location (j,k).
6) It is unallowable for the extent of groundwater contamination to increase beyond initial conditions at any time during the remediation.
7) Total pumping and infiltration rates must be balanced at the beginning of every management period
where:
is the number of injection wells operating in year i.
In summary, the optimization problem can be stated as follows: find the combination of simulated extraction and injection well locations and their operating rates that minimize the cost of reducing RDX and TNT concentrations within a 20 year time horizon while satisfying all the constraints on well, remedy operations, and plume behavior. The remedial system designs use the calendar year 2003 as the starting point.
The groundwater flow and transport models provided for the study simulate 20 years with four management periods of five years each. The extraction and injection flow rates can vary across the management periods, but not within a management period. The locations of the extraction and injection infrastructure can only vary between candidate solutions, not between management periods.
Optimization solution approach
The optimization problem is solved by defining the regions to search for the globally optimal settings of the decision variables, and prescribe how to conduct the search. The decision variable search regions mimicked, as closely as possible, the regions established for new well locations and infiltration basins by the MGO team. The MGO team confined their search for new extraction well locations to regions where the plume densities were highest as shown on Figure 2. Each region contains the areas of the highest concentration for the three principal lobes of the two contaminant plumes. In addition to these search areas for extraction wells, locations for three new infiltration basins, IFA, IFB, and IFC in the original study at the extent of the RDX plume to the east, southeast, and southwest as alternatives to the four preexisting recharge basins. In all, a total of 11 candidate regions exist for locating decision variables: four extraction areas and seven injection/recharge areas; the locations as defined by the MGO team. These 11 regions make up the infrastructure search areas of the remedy used to alter the distribution of RDX and TNT in groundwater. Active remediation solutions require nonzero total extraction and injection fluxes; hence, there must be at least one extraction location and one injection location active for all viable candidate solutions.
Within the extraction locations, a well(s) can be located anywhere in the search box (we restrict the well position to a model node). Table 2 provides the box location and number of candidate position locations for one or two wells in each singular box.
Table 2. Number of candidate extraction well locations in each of the three search boxes
Allowing an extraction well(s) to be located in any of the three search boxes results in 302,253 potential extraction well location combinations. The 128 different candidate location configurations for the infiltration basins results in a total of 38,688,384 candidate infrastructure system designs. Simultaneously, when determining the infrastructure configuration design, the water flow rates are optimized. This magnitude of options illustrates the difficulty of finding the optimal solution either by random searching or by using the SME subjective engineering judgment.
The approach used to solve the Umatilla problem consisted of the following generalized automated search strategy:
1) Evaluate the cost of the RIP. The RIP consists of three extraction wells and three infiltration basins. Store these results as a current minimum.
2) Begin evaluation of different combinations of extraction wells and infiltration basins.
a. Set number of evaluation epochs equal to one.
a. Initial extraction location: Initiate the solution using the study maximum number of new wells (two), located in Box 1 (which contains the RDX and TNT plumes).
a. Initial injection location: IFL located in Box 1, the innermost infiltration basin.
a. Initiate search strategy. Begin to cycle through and test the 4 pumping areas and 7 injection areas for quality of solution, use extraction rate total = system maximum capacity of 1170 (gpm) (which is 1,300 (gpm) * 90% uptime) as done by the MGO team.
a. Set search cycle a minimum number of evaluations (12).
a. Upon finding a cost lower than current minimum cost:
b. Select the solutions with the two lowest costs.
b. Explore these regions by conducting brief local random searches (LRS) on each solution to (statistically) determine the most promising configuration of wells and basins. Perform initial analysis using the physicallybased model. The optimization search progresses and generates an increasing number of function evaluations. Machine learning is invoked to provide candidate solution objective function evaluation estimation which reduces computation burden when model evaluations models are extensively time consuming.
a. Store these results (configuration, rates, cost).
a. Evaluate best current solution.
b. If solution improves, replace previous remedial design configuration and objective function value.
b. Else, continue.
3) Initiate a global optimization analysis using the LGO search options of GARS followed by GRG; initialize using currently best identified solution.
a. Upon finding a cost lower than current minimum cost:
i. Store these results.
a. Else, continue.
4) If termination criteria met,
a. Stop.
a. Else, increment number of evaluation epochs by 1 until the userspecified maximum number reached. Store and print results.
a. Go to step 2.
PBMO’s partition and exploration approach enables implementation of the search strategy be conducted on multiple central processing units (CPUs). However, this test used a sequential implementation to mimic ESTCP test conditions. The search continues until the termination criterion satisfied. The termination criterion can either be the targeted optimal cost, the total number of flow and transport simulations, the number of simulations since the optimal value was last improved, or the total simulation (CPU) time consumed. In this case, the termination criterion was the optimal total remedial cost as published by the MGO team. The initial starting point for the total flow rate is the maximum treatment plant flow rate 1,170 (gpm) and mimics the MGO team. Initially, all extraction wells in Box 1 equally pumped; all extracted water injected into the single infiltration gallery IFL, and new candidate extraction wells placed randomly during the infrastructure search phases.
Results and discussion
The total net present value cost from PBMO ($1,664,085) significantly improved upon the costs committed for remediation at the site by implementing the RIP ($3,836,285) and the design achievable via SME ($2,230,905). The PBMO results mimicked MGO ($1,664,395) and SOMOS ($1,663,841). Table 3 presents the well location strategies and cost results for RIP, SME, MGO and PBMO (since MGO and SOMOS were similar).
Table 3. Optimal pumping strategies found using SME, MGO and PBMO for formulation 1 at Umatilla compared with existing RIP design
The termination criterion in the test was the MGO cost value. PBMO found a lower cost that MGO in fewer than 120 model evaluations. Additional optimization analysis performed during reliability testing produced multiple solutions with lower costs  as low as $1,663,240 – and with different new well extraction locations and different rates in the extraction wells. These subsequent findings illustrate the multiextremal structure of this design problem, thereby calling for global optimization based solution approaches.
Simulation model reliability over the range of feasible inputs is essential for determining the globally optimal variable values. We observed that 10.7% of the viable candidate designs simulated in the groundwater flow and transport models during the optimization search failed to converge. This primarily occurred when the extraction well flow rates were set near the high end of the acceptable range and which caused flow model cells to dewater. In instances where model cells become dewatered, the groundwater flow simulator could not converge to a solution without the application of a nonphysical lower bound on the head in the dewatered cell^{a}. Given that the two simulators execute sequentially during a function evaluation (groundwater flow followed by fate and transport), the failure of the flow simulator to provide a convergent solution prevents the transport simulator from predicting contaminant distributions resulting in a lost candidate design evaluation cycle.
PBMO addresses issues that arise in model simulations due to these harsh modeling scenarios imposed by formal optimization. PBMO examines the model simulation solution time and the flow and transport mass balance errors. Automated solver parameter adjustments can take place if the solution becomes unstable or inefficient. If a model nevertheless fails to converge after a userspecified maximum number of numerical solver parameter setting attempts, penalty functions divert the optimal search from exploring this solution region. Hence, while the algorithm handles the nonconvergent model issue, should the simulations be noted to fail to converge either at an appreciable rate or in regions of the search space where the suspected location of the optimal value, a more robust model code should be used. This will alleviate the risk of not locating the true globally optimal solution. This study design necessarily used the modeling system used by the other teams in spite of the 10.7% model simulation failure rate.
Table 4 presents a breakdown of the capital, and operations and maintenance costs for each of the four strategies. Figure 4 shows the well locations for RIP and T&E conducted by the SME’s, and the optimal well locations determined by MGO and PBMO The PBMO solution places two new wells into the TNT plume, as did MGO.
Table 4. Breakdown of the capital and O&M costs of RIP, SME, MGO and PBMO designs
Figure 4. Well location configurations generated by various techniques: RIP, PBMO, MGO and SME (using trial and error).
Examining the performance comparisons of the algorithms^{b} reported in the ESTCP study for the number of flow and transport simulations, PBMO found its solution using under 120 flow and transport simulations as compared to the estimated 5,000 simulations reported by the MGO investigators [see Volume II of the DoD/ESTCP report (Minsker et al., 2004)]. This is a significant improvement in efficiency, by an estimated factor of over 40. Furthermore, PBMO ran unattended, whilst MGO required numerous human interventions.
The authors believe that the observed increased efficiently lies in the integrated optimization search strategy. The PBMO approach supports rapid determination of “good” solution spaces that work synergistically with LGO, the core global optimization algorithm. The investigation of (Rios and Sahinidis 2012) supports this position via the results of testing 23 optimization algorithms, including LGO. The results of that comparative study indicate LGO to be orders of magnitude more efficient and reliable than the techniques selected by the original ESTCP study teams (see Figure 5, adapted from that study). Regarding the final solution to the Umatilla problem, the solutions generated by the optimization techniques are conceptually similar: all three utilize the same existing extraction wells, EW1 and EW3; both use the same infiltration basins, IF2 and IF3; and both locate new extraction wells near to the center of mass of the TNT plume, and remediate the site in 4 years. Figures 6, 7, 8, 9 show the RDX and TNT plumes at the end of Year 1 through Year 4― the time of remediation completion. Figure 10 and Figure 11 illustrate the maximum concentration of both contaminants over time in the top layer of the model. The (slight) over design by MGO compared with PBMO is evidenced regarding the TNT remediation results. The TNT remediation did not occur at the same time as the RDX remediation and the water quality was remediated cleaner than required by the project specifications.
Figure 5. Comparison of selected optimization techniques.
Figure 6. Extent of RDX and TNT for PBMO and MGO optimal solutions after year 1.
Figure 7. Extent of RDX and MGO optimal solutions after year 2.
Figure 8. Extent of RDX and MGO optimal solutions after year 3.
Figure 9. Extent of RDX and MGO optimal solutions after year 4.
Figure 10. Calculated maximum concentrations of RDX in the shallow aquifer model layer 1 (Starting at the end of 2002).
Figure 11. Calculated maximum concentrations of TNT in the shallow aquifer model layer 1 (Starting at the end of 2002).
Remediation well starting position: sensitivity tests
The algorithm’s reliability test consisted of assessing its ability to converge on the optimal result from different illplaced starting locations. The well location search region used was the same search region defined for the northernmost box, the one that contains the initial RDX and TNT plumes as shown in Figure 2. However, instead of using the best configuration of the extraction well locations determined thus far (i.e. from the infrastructure search), the search box corners specify the initial new extraction well locations. The test scenarios consisted of six combinations of initial locations for the two new extraction wells at the four corners. Figure 12 shows the six different initial starting configurations for the two new extraction wells represented as the green triangles on the various corners. The results of interest are the final solution, the optimal placement of two new extraction wells to accompany EW1 and EW3, the number of simulation/optimization iterations required, and the elapsed time needed to achieve the solution. Table 5 presents the test results. Each of the six test cases found the same optimal locations. The number of model simulations varies between 124 and 127. In every scenario, identification of the optimal solution occurred in less than 3 hours of CPU time. In nearly 100 iterations, the GARS optimization algorithm found the optimal solution regardless that the search initiates from locations outside the contaminant plume. Reliability and minimal dependence on the starting solution initialization are fundamental and desirable features of a highquality global solver.
Figure 12. Starting locations of new candidate wells for reliability testing of PBMO.
Table 5. PBMO computational performance during robustness testing
Conclusions
The study demonstrates the effectiveness of an automated, costeffective, and broadly applicable approach to groundwater remedial design. The application of the PBMO methodology to the Umatilla case study site represents a rigorous testing and validation exercise. This numerical analysis efficiently and automatically identified the best solution found by others. Extension of the approach for evaluation and optimal aquifer remediation management for different sites in different geologic settings is accomplished by using a site specific calibrated groundwater flow and transport model. The approach identified the optimal solution about 40 times faster than any of the other methods by reducing the number of time consuming flow and transport model evaluations. Comprehensive automation promotes efficiency, effectiveness and usability. Merit for reliably optimizing complicated systems is achieved through the ability to overcome the inherent nonconvergence difficulties that occur when using numerical models for process simulation.
Endnotes
^{a}This approach was explicitly applied to one of the other study sites in the DoD/ESTCP study – see Section 3.2.3.4 in Volume I of Minsker et al. (2004).
^{b}The model files representing the optimal solution identified by MGO and SOMOS are on the web site, however, the input files to MGO or SOMOS to recreate the optimization study and generate the optimal solution are not available. Therefore, it was not possible to independently verify the reported operational performance of MGO or SOMOS. Hence, we use the reported values for the number of model runs, with MGO being the lesser number of ~5,000. To put this efficiency result in perspective, it would take nearly 5 days of clock time on a modern computer (~1.4 minutes per flow and transport simulation) to solve the problem using MGO, whereas PBMO found a somewhat better solution in less than 3 hours. In addition, the MGO solution required several stops and starts as well as other interceding actions by the investigators to reach the final answer presented in the study. PBMO required a single execution without any human intervention: the optimization process was completely automated.
Abbreviations
BB: Branch and Bound (optimization method in LGO); CCW: Capital Costs of New Wells; CCB: Capital Costs of New Infiltration Basins; CPU: central processing unit; DoD: U.S. Department of Defense; EPA: U.S. Environmental Protection Agency; ESTCP: Environmental Security Technology Certification Program; EW: extraction well; FCL: Fixed Costs of Labor; FCE: Fixed Costs of Electricity; GA: Genetic algorithm; GAC: Granular activated carbon; GARS: Globally Adaptive Random Search (optimization method in LGO); gpm: Gallons per minute; GRG: Generalized reduced gradient (optimization method in LGO); HGL: HydroGeoLogic, Inc.; kg: Kilograms; LGO: Lipschitz Global Optimization software package; LP: Linear programming; LRS: Local random search; μg/L: Micrograms per liter; MGO: Modular Groundwater Optimization; MODFLOW: Modular FiniteDifference GroundWater Flow ModelMS, multi start random search (optimization method in LGO); MT3DMS: Modular 3D MultiSpecies Transport Model; NPL: National Priorities List; NPV: Net present value; OA: Outer approximation; O&M: Operations and maintenance; PBMO: PhysicsBased Management Optimization; ppb: Part per billion; RAM: Randomaccess memory; RDX: Royal Demolition Explosive; RIP: Remedy In Place; SA: Simulated annealing; SLA: Sequential linear approximation; SLP: Sequential linear programming; SME: Subject matter expert; SOMOS: Simulation/Optimization Modeling Optimization Software System; SQP: Sequential quadratic programming; TNT: 2,4,6Trinitrotoluene; TS: Tabu search; Umatilla: Umatilla Chemical Depot; VCE: Variable Costs of Electricity for Operating Wells; VCG: Variable Costs of Changing GAC Units; VCS: Variable Costs of Sampling.
Competing interests
The authors declare that they have no competing interest.
Authors’ contributions
The work presented here was carried out in collaboration between all authors. LD defined the research theme, designed and implemented the overall optimization approach. JP provided LGO optimization component and TL coded the objective function calculator. All authors have contributed to the paper preparation, and have seen and approved the manuscript.
Authors’ information
Larry M Deschaine PE is a graduate of MIT (1984), the University of Connecticut (1992) and is completing his PhD work in Complex Systems at Chalmers University of Technology in Göteborg, Sweden (expected 2013). He is a Principal Engineer at HydroGeoLogic and a recognized expert in the development and application of simulation and optimization techniques to largescale energy and environmental engineering problems. His work covers the range of surface water, ocean and groundwater characterization and simulation; remedial design and construction; monitoring optimization; detection of unexploded ordnance, and; the development of selflearning adaptive algorithms which support optimization of electric power distribution networks including grid integration of renewables. His work has received numerous awards, including a US VicePresidential Hammer Award, and he has written more than 100 journal articles, book chapters, and other technical documents.
János D. Pintér is a researcher and practitioner with four decades of experience in the area of systems modeling and optimization, including algorithm and software development. He holds MSc, PhD and DSc degrees in Mathematics, with specialization in the Operations Research area. He has authored and edited books including an awardwinning monograph, and has written more than 200 journal articles, book chapters, and other technical documents. Dr. Pintér is an active member and officer of several professional organizations, and he serves on the editorial board of several professional journals.
Theodore P. Lillys PE is a Senior Engineer at HydroGeoLogic, Inc. with over fifteen years of experience in research and applications in groundwater modeling, optimization, statistics, and software development. He holds and MS (1997) in civil and environmental engineering specializing in subsurface hydrology, numerical methods, and mathematical optimization techniques. Mr. Lillys’ areas of expertise lie in groundwater modelling, management optimization of subsurface remedial designs and resource management, optimization of long term monitoring networks, and software development.
Acknowledgements
This research was conducted in concert with LD’s PhD studies in Complex Systems at Chalmers University, with financial support by HydroGeoLogic, Inc, 11107 Sunset Hills Road, Suite 400, Reston, VA 20190 (http://www.hgl.com). An animation of this case study has been prepared; see (http://www.hglsoftware.com/cleanup.cfm).
References

Ahlfeld DP, Barlow PM, Mulligan AE (2005) GWM – a groundwater management process for the US Geological Survey modular groundwater model (MODFLOW – 2000). US Geological Survey.
OpenFile Report 20051072, 124 Pages, 2005. https://water.usgs.gov/nrp/gwsoftware/mf2005_gwm/OFR2005_1072.pdf webcite
PubMed Abstract  Publisher Full Text 
Deschaine LM (1992) Cost evaluation and optimization of ground water pump and treat programs, MSCE thesis. Storrs, Connecticut: University of Connecticut.

Deschaine LM, Ahlfeld DP, Ades MJ, O’Brien D (1998) An optimization algorithm to minimize the life cycle cost of implementing an aquifer remediation project  theory and case history. Society for Modeling and Simulation International, Simulators International XV, Simulation Series 30(3):5358

Deschaine LM, Regmi S, Patel JJ, Fox TA, Ades MJ, Katyal A (2001) Design optimization of groundwater quality management challenges using the outer approximation method, The Society for Modeling and Simulation International. Seattle, WA, USA: Advanced Simulation Technology Conference. pp pages 8893
April 2001

Deschaine LM (2003) Simulation and optimization of largescale subsurface environmental impacts; investigations, remedial design and long term monitoring. Journal of Mathematical Machines and Systems, Kiev 3(4):Pages 201218

Deschaine LM, Pintér JD (2003) A comparison of the outer approximation method and lipschitz global optimization on optimal groundwater quality management. Atlanta, Georgia: Presented at the INFORMS Annual Meeting Conference.
October 1922, 2003

Deschaine LM, Nordin JP, Pintér JD (2011) A computational geometric / information theoretic method to invert physicsbased mecmodel attributes for mec discrimination. Journal of Mathematical Machines and Systems, National Academy of Sciences of Ukraine, Kiev 2:5061

Guo X, Zhang CM, Borthwick JC (2007) Technical Report. Water Resources Research Journal 43(No. 8):WO841614 Pages
August. http://onlinelibrary.wiley.com/doi/10.1029/2006WR004947/abstract webcite

Harbaugh AW, McDonald MG (1996) User’s documentation for modflow96, an update to the u.s. geological survey modular finitedifference groundwater flow model. Reston VA: U.S. Geological Survey OpenFile Report 96–485. pp 56485 PubMed Abstract

HydroGeoLogic, Inc. (HGL) (2012) MODFLOWSURFACT™ Version 3.0; User’s manual and guide. VA, USA.
20190
PubMed Abstract 
ITRC (2007) InSitu Bioremediation of Chlorinated Ethene DNAPL Source Zones: Case Studies. Prepared by The Interstate Technology and Regulatory Council, Bioremediation of Dense NonAqueous Phase Liquids (Bio DNAPL) Team. Refer to Chapter 9 “Simulation and Optimization of Subsurface Environmental Impacts; Investigations, Remedial Design and Long Term Monitoring of BioNAPL Remediation Systems”.. 128147
April 2007. http://www.itrcweb.org/Guidance/GetDocument?documentID=11 webcite

McDonald MG, Harbaugh AW (1988) A modular threedimensional finitedifference groundwater flow model, Techniques of Water Resources Investigations Book 6. Washington, D.C: U.S. Geological Survey.

Minsker B, Zhang Y, Greenwald R, Peralta R, Zheng C, Harre K, Becker D, Yeh L, Yager K (2004) Final technical report for application of flow and transport optimization codes to ground water pumping and treat systems, Technical Report to the Environmental Security Technology Certification Program. Volumes IIII. TR2237ENV. California: Engineering Service Center, Port Hueneme.
http://www.frtr.gov/estcp/estcp.htm webcite

Peralta RC (2004) “SOMOS: Simulation/Optimization Modeling System”. User’s Manual, Software Engineering Division, Department of Biological and Irrigation Engineering. Logan, UT: Utah State University. p 48 Pages
April 2004. http://www.frtr.gov/estcp/source_codes.htm webcite

Pintér JD (1996) Global optimization in action. Kluwer Academic Publishers, Dordrecht Boston London, 1996. New York: Now distributed by Springer Science + Business Media.

Pintér JD (2002) Global optimization: software, test problems, and applications. In: Pardalos PM, Romeijn HE (eds) Handbook of Global Optimization, Volume 2, Dordrecht: Kluwer Academic Publishers. pp 515569

Pintér JD (2009)Pardalos PM, Coleman TF (eds) Software development for global optimization, Providence, RI: American Mathematical Society.

Pintér JD (2013) LGO  A model development and solver system for globallocal nonlinear optimization user’s guide. Canada: Published and distributed by Pintér Consulting Services, Inc.
http://www.pinterconsulting.com webcite. First edition: June 1995; Current edition: March 2013

Rios M, Sahinidis N (2012) “Derivativefree optimization: A review and comparison of software implementations”. J Glob Optim.
http://link.springer.com/content/pdf/10.1007%2Fs108980129951y.pdf webcite
Publisher Full Text 
Systems Simulation/Optimization Laboratory (SSOL) (2002) Simulation/Optimization Modeling System (SOMOS) Users Manual. Logan, UT: SS/OL, Biological & Irrig. Eng. USU. p 457

U.S. Army Corps of Engineers (USACE) (1994) Defense Environmental Restoration Program, Final Record of Decision, Umatilla Depot Activity Explosives Washout Lagoons Ground Water Operable Unit.
June 7, 1994

U.S. Army Corps of Engineers (USACE) (1996) Final Remedial Design Submittal, Contaminated Groundwater Remediation, Explosives Washout Lagoons, Umatilla Depot Activity, Hermiston Oregon.
January 1996

USEPA (2012) National strategy to expand superfund optimization practices from site assessment to site completion.
(September 2012) (PDF) OSWER 9200.375 http://www.epa.gov/oerrpage/superfund/cleanup/postconstruction/optimize.htm webcite

Zheng C, Wang PP (1999) MT3DMS: A modular threedimensional multispecies transport model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems; documentation and user’s guide, Contract Report SERDP991. Vicksburg, MS: U.S. Army Engineer Research and Development Center.
available at http://hydro.geo.ua.edu/mt3d webcite

Zheng C, Wang PP (2002) MGO – A modular groundwater optimizer incorporating modflow/mt3dms, documentation and user’s guide. Draft.
April 2002. http://www.frtr.gov/estcp/source_codes.htm webcite
PubMed Abstract 
Zheng C, Wang P (2003) “Application of flow and transport optimization codes to groundwater pumpandtreat systems: Umatilla Army Depot, OR”. Technical Report, Revised version 2/2003. Tuscaloosa, AL: University of Alabama. p pp. 41
http://www.frtr.gov/estcp/demonstration_sites.htm webcite