Hydraulic gradient control for groundwater contaminant removal

Journal of Hydrology, 76 {1985) 85--106 Elsevier Science Publishers B.V., Amsterdam -- Printed in The Netherlands

85

[2] HYDRAULIC GRADIENT CONTROL FOR GROUNDWATER CONTAMINANT REMOVAL

DOROTHY FISHER ATWOOD and STEVEN M. GORELICK

U.S. Geological Survey, Menlo Park, CA 94025 (U.S.A.) (Received April 24, 1984; accepted for publication June 13, 1984)

ABSTRACT Fisher Atwood, D. and Gorelick, S.M., 1985. Hydraulic gradient control for groundwater contaminant removal. J. Hydrol., 76: 85--106. The Rocky Mountain Arsenal near Denver, Colorado, U.S.A., is used as a realistic setting for a hypothetical test of a procedure that plans the hydraulic stabilization and removal of a groundwater contaminant plume. A two-stage planning procedure successfully selects the best wells and their optimal pumping/recharge schedules to contain the plume while a well or system of wells within the plume removes the contaminated water. In stage I, a combined groundwater flow and solute transport model is used to simulate contaminant removal under an assumed velocity field. The result is the approximated plume boundary location as a function of time. In stage II, a linear program, which includes a groundwater flow model as part of the set of constraints, determines the optimal well selection and their optimal pumping/recharge schedules by minimizing total pumping and recharge. The simulation--management model eliminates wells far from the plume perimeter and activates wells near the perimeter as the plume decreases in size. This successfully stablizes the hydraulic gradient during aquifer cleanup.

1. INTRODUCTION In r e c e n t y e a r s t h e i n c r e a s e in t h e n u m b e r o f i n c i d e n t s o f g r o u n d w a t e r p o l l u t i o n h a s l e a d t o a g r e a t e r i n t e r e s t in t h e b e s t w a y t o m a n a g e a n d r e s t o r e c o n t a m i n a t e d areas. A c c o m p a n y i n g t h i s i n t e r e s t has b e e n t h e d e v e l o p m e n t of aquifer management tools which integrate groundwater flow and solute t r a n s p o r t s i m u l a t i o n w i t h o p t i m i z a t i o n m e t h o d s . As d e m o n s t r a t e d here, s u c h m e t h o d s c a n aid t h e h y d r o l o g i s t in e v a l u a t i n g a n d p l a n n i n g s c h e m e s f o r a q u i f e r r e s t o r a t i o n . In p a r t i c u l a r , we p r e s e n t a d e s i g n m e t h o d o l o g y f o r h y d r a u l i c g r a d i e n t c o n t r o l a i m e d at c o n t a i n i n g a n d r e m o v i n g g r o u n d w a t e r contaminants. A l t h o u g h m a n y works have c o m b i n e d aquifer simulation with optimiz a t i o n t e c h n i q u e s s u c h as l i n e a r a n d q u a d r a t i c p r o g r a m m i n g , f e w s t u d i e s have considered t h e p r o b l e m o f c o n t a m i n a n t p l u m e c o n t a i n m e n t . Molz and Bell ( 1 9 7 7 ) c o n s i d e r e d t h e p r o b l e m o f h y d r a u l i c s c r e e n i n g t o k e e p a z o n e o f

86

artificially emplaced hot groundwater from migrating. They applied the embedding method of Aguado and Remson (1974). In that method the groundwater flow equation is included in a linear programming (management) model as constraints along with other limitations on pumping rates and local hydraulic gradients. The application successfully identified pumping and recharge rates from wells at fixed locations. The limitation of the method is that in practice it can only be applied to small steady-state flow problems. Both excessive computer storage requirements and numerical difficulties make the method infeasible for field-scale transient problems. [See the review paper by Gorelick (1983) for a discussion of this point.] Remson and Gorelick (1980) demonstrated the use of the embedding approach to identify well locations and rates to keep a contaminant plume from migrating. Other restrictions included local dewatering requirements and export water limitations. Another powerful technique in groundwater management modeling is the use of a response matrix, generated from a groundwater flow model, in conjunction with linear or quadratic programming. The matrix shows the influence of pumping or recharge wells at potential well sites upon drawdowns at specified observation locations. This approach has been applied to a variety of hydraulic management problems involving maximization of well yields and minimization of production costs (Lee and Aronofsky, 1958; Deininger, 1970; Wattenbarger, 1970; Maddock, 1972; Rosenwald and Green, 1974}. A management model based on a response matrix has yet to be applied to problems involving groundwater contaminant control. We describe a procedure that integrates the use of linear programming, groundwater flow simulation and solute transport simulation to solve the latter problem.

2. P R O B L E M D E S C R I P T I O N

A procedure is demonstrated here that selects the best wells and their pumping/recharge schedules to contain a contaminant plume, while a well or system of wells within the plume removes contaminated water. Wells used for hydraulic gradient control are picked from preselected well sites and must be outside of the zone of contaminated groundwater. Pumping at upgradient wells or recharging at downgradient wells produces a flat or inward gradient at the plume boundary. As the plume decreases in size, potential gradient control wells initially within the zone of contamination (and therefore inactive} become surrounded by clean water. Once the wells are outside the contaminant plume they may be selected by the management model to help control the hydraulic gradient. Likewise, well sites distant from the plume boundary in later periods may be shut down in the management solution. The hydraulic gradient control design procedure (Fig. 1} is divided into two stages to linearize the water quality management problem. In stage I, solute transport simulation is used to predict the location of the shrinking

87

GROUNDWATER CONTAMINATION PROBLEM Assumed Velocity Field

STAGE I L TRACK PLUME BOUNDARY UNDER NEAR ZERO GRADIENT I

"x#

/

Approximated Contamination Distribution {D

CD Z

STAGE II

IDENTIFY OPTIMAL WELL SELECTION AND RATES

tO {D m

Known Velocity Field

CHECK SOLUTION

L

Fig. 1. H y d r a u l i c g r a d i e n t c o n t r o l p r o c e d u r e .

plume boundary over time, assuming that the regional hydraulic gradient has been effectively flattened in the vicinity of the contaminant plume {sometimes called hydraulic screening). The geometry of the contaminant plume, determined in stage I, is used to identify potential well locations. The optimal well selection and pumping/recharge schedules are given b y a simulation--management model in stage II. Finally, the resulting optimal schedules of stress rates are used in a simulation model o f groundwater flow and solute transport to verify the solution and give a detailed view of contaminant removal and redistribution over space and time.

3. M O D E L A R E A

The R o c k y Mountain Arsenal near Denver, Colorado, U.S.A., is chosen as the setting for this study because o f its well-known geology and hydrology

88 MODEL

AREA

£"~" ~ / i ~,'#

........ ~

I

i ARSENAL i

if,--,_.

,~I

i!

i"" ROCKY ! MOUNTAIN I

IF

,~- . . . . .

-i.._J

r-- .r-

(DENVER.I~~

! x

r'4 .

\

"

h

j

L,I_

'"- ~

0

5

10

I

I

I

15 Kilometers

COLORADO

I

Fig. 2. Location map.

(Fig. 2). The Arsenal is a military facility that was designed to manufacture and process toxic chemicals for military and industrial use. Poor disposal techniques have resulted in a long history of groundwater contamination. The specific area considered here straddles the north boundary of the R o c k y Mountain Arsenal where preventing contaminant migration has been crucial. Flow in the unconfined aquifer is generally from south to north. The goal of the restoration procedure is to prevent further contaminant migration during aquifer cleanup, especially b e y o n d the northern boundary o f the Arsenal. The model area is a small part (3 km 2 ) of a large groundwater basin. This was discretized into 702 finite-difference nodes (445 active) each 76.2 m on a side {Fig. 3). The boundaries were approximated. On the southwest side a flux boundary was constructed based upon known hydrologic parameters and flows (Konikow, 1977). This was done by running the steady-state simulation with constant-head boundaries and then calculating the resulting fluxes across the southwest boundary. The northern constant-head boundary is based on steady-state water levels. Most of the west boundary and a small portion in the southeast are areas o f high bedrock where the alluvium is usually unsaturated. These are represented by no-flow boundaries. The eastern no-flow boundary is a flow line which divides the flow to the east from flow in this area.

89

- ~ I.Z 1~" I x

...... ~" N o - F l o w Boundary === =-- Constant-Head Boundary = = o - - , Flux Boundary

~

Contaminant Plume Central Contaminant Removal Well O Well Locations . ~ o 7 0 " Water Table Elevation (Meters above Sea Level) •

Fig. 3. Model area. Previous studies (Konikow, 1975, 1977) have used subsurface data in and to the north o f the Arsenal to describe the hydrogeologic characteristics of the alluvial aquifer. The hydrogeologic parameters determined in these studies are the basis for the groundwater and solute transport models used in this study. The steady-state hydraulic head distribution (Fig. 3) generally matches the natural water levels presented in the studies of Konikow (1975, 1977). There is a shallow gradient in the south which steepens to the north. The transmissivity is extremely variable (Fig. 4). In general, the transmissivity is higher in the southeastern part of the model area. Saturated thickness

90 ?-

r---T--

,L r

.

.

.

.

\-

/; "~-

ti.-

/

0

/

500 Meters

~

//1!.",..,

/

~

I

I l l l l l '~1"~ ,';" N o - F l o w Boundary

f ~5 ~ Transmissivity (in m2/day) Fig. 4. Transmissivity distribution.

varies from 0 to 1 0 m . Constant values of effective porosity (0.30) and storage coefficient (0.25) are based on the estimates by Warner (1979). The longitudinal and transverse dispersivities are each 30 m as determined b y Konikow (1977). Groundwater contamination is an existing problem at the R o c k y Mountain Arsenal. F o r purposes of exposition in this study, however, a hypothetical non-reacting contaminant was introduced into the model area using solute transport simulation. The boundary o f the plume is defined at the 250-mg-1-1 contour which is taken as the "safe" concentration. 4. GOVERNING EQUATIONS

The hydraulic gradient control procedure requires the use of models of groundwater flow and solute transport. Groundwater flow is simulated using

91

the finite-difference model developed by Trescott et al. (1976). The equation for a confined aquifer is used as an approximation of the unconfined aquifer at R o c k y Mountain Arsenal. As long as the drawdowns are small relative to the saturated thickness this approximation will be fairly good. The governing equation for two-dimensional flow in a heterogeneous isotropic aquifer is:

a

axi

T(x) 0h

= S-- +W, at

i,j=1,2

(1)

where T(x) = transmissivity [L 2 T - I ] ; h = hydraulic head [L] ; S = storage coefficient [dimensionless] ; t ~: time [T] ; W = volume flux per unit area [L T -1 ] (source/sink term); and xi,xj = space variables [L]. Solute transport is combined with groundwater flow simulation in the c o m p u t e r code o f Konikow and Bredehoeft (1978). Once the groundwater velocity field is determined using Darcy's law and the dispersion coefficients are determined using the identities of Scheidegger (1961), the method-ofcharacteristics is used to simulate solute transport based on the advection-dispersion equation (transport equation):

a(Cb) at

-

3( a_~__xi) axi bDij

a(bCVi) axi

C'W e

,

i, j = 1 , 2

(2)

where C = concentration [M L - 3 ] ; Vi(W) = velocity [ L T - 1 ] ; DiI(V) = coefficient o f hydrodynamic dispersion [L: T - l ] ; W = volume flux per unit area [L T -1] (source/sink term); C ' = concentration of source/sink fluid [M L -3 ] ; b = saturated thickness [L] ; e = porosity [dimensionless] ;xi, x1 = space variables ILl ; and t = time IT]. Simulation of contaminant transport and well field optimization for hydraulic gradient control cannot be combined into a single linear management problem because o f system nonlinearities. This is because the pumping/recharge rates, W, are decision variables determined in the management problem. With these rates initially unknown, the groundwater velocity field, V, is also u n k n o w n and the simulation of contaminant transport is impossible. These nonlinearities result from multiplying u n k n o w n velocities and concentrations as shown in the first and second terms on the right-hand side of the governing advection--dispersion equation (2). Our two-stage procedure deals effectively with the above nonlinearities.

5. H Y D R A U L I C G R A D I E N T C O N T R O L D E S I G N P R O C E D U R E

5.1. Stage I: Simulate contaminant distribution In stage I contaminant transport is simulated through time using an a~l,med velocity field (corresponding to successful hydraulic screening)

92

to determine the location of the plume boundary relative to the potential gradient control wells. Three steps are necessary in stage I: (1) an approximate velocity field is assumed; (2) the central contaminant removal system is designed; and (3) the contaminant distribution (removal) is simulated. Since the true velocities will remain u n k n o w n until the management problem has been solved, the first task in stage I is to determine an appropriate velocity field approximation. Such a velocity field can be determined by noting the effects of successful hydraulic gradient control. In the managed system, wells outside the plume perimeter maintain a zero or inward gradient while the contaminant removal system (pumping wells) operates. Therefore, a good approximation of the velocity field is produced b y an initially flat hydraulic surface, affected only b y the contaminant removal system. ~

n a ~ m m a l i i m u u m o immimmlmanln

inimau

~

~

u

n

~

muammmJ --~ /

/~

;/ \-,

~

5

0

\,

}

N

0 I I '"(:"-'";; ... m oB ~ ~Lc~O-

500 Meters ~ _

No-Flow Boundary Equal Value Constant-Head Boundary Trial Well Location Concentration in mg/£

Fig. 5. Initial c o n t a m i n a n t c o n c e n t r a t i o n d i s t r i b u t i o n w i t h trial well locations.

93 The speed at which the contaminant is reduced to a safe level (defined as 2 5 0 m g l -]) depends on the number, the placement, and the pumping rates o f the contaminant removal wells. In the illustrative problem, we use a single decontamination well pumping at the fixed rate of 2.83 1 s -]. The most effective location was determined by trial and error with simulations using four different well locations (Fig. 5). The results of the simulations are shown in Fig. 6. Well A is located at the node with the highest initial contaminant concentration. This was considered a good initial location for the contaminant removal well. Cleanup using this well, however, requires almost twenty-two years, the longest time of all the trial locations. Studying the contaminant distribution shows that the contaminant is more effectively removed from the southern portion of the plume. High concentrations remain in the north much longer than the south. To counteract this effect, the other trial wells were located north o f well A. When pumping from well C, located 154.2 m north of well A, the concentration falls below 250 mg 1-1 in the entire system in less than 13 yr. Pumping in the other two locations produces intermediate results. Well C was chosen to be the location for the contaminant removal well because it restores water quality most effectively. 1750

[

I

I

I

kk\

A

1250 E z o I< n-. ~z IJJ o z o £3

750

250

- - Safe Level 12"WSe~'leCrs l~irSsYe&r D

21"5elYe;rs/

0 0

5

10

15

20

25

YEARS

Fig. 6. Maximum contaminant concentration in plume when pumping at indicated trial locations.

94

Redistribution o f the contaminant through time using well C is shown in Fig. 7. Constant-head boundaries, equal in value, replace the flux boundaries creating a system with no initial hydraulic gradient. All the other parameter distributions remain the same, e.g., transmissivity, saturated thickness and porosity. The center well pumps at a constant rate of 2 . 8 3 1 s -1, creating an inward gradient. For the first 10yr., the contaminant is removed more quickly from the southern portion of the plume as shown by the wider spacing between the 0.5-yr. contours. The entire plume is removed in just over 1 2 y r . This approximation of the plume geometry is used in the design of the R o c k y Mountain Arsenal management problem in stage II.

~ I+

II0 ':("-""~; "---• j8 ~-

m

m

m

m

m

l

l

m

m

m

m

m

m

m

Y

% • tx \-,

"<% 500 Meters ~'~

-

No-Flow Boundary Equal Value Constant-Head Boundary Central Contaminant Removal Well Plume Boundary Contours by Year

Fig. 7. Simulated plume boundary with time.

95 5. 2. Stage H: Select best wells and rates Linear programming combined with groundwater flow simulation is used to select the wells and determine the pumping/recharge schedules that most effectively stop the migration of the plume b y controlling the hydraulic gradient during cleanup. This is accomplished with a single global optimization which minimizes the total pumping/recharge over all time subject to a series o f hydraulic constraints that force the gradient to be toward the center contaminant removal well. Objective: Min {e}W{N} q- {e}T{v}

(3)

Constraints: [R]{u} -- [R]{v} < {g}

(4a)

uc = 2.83 l s -1

(4b)

with

where [R] = matrix of gradient response coefficients at gradient check locations; { u} = pumpage vector; ue = pumping rate at central removal well; { v} = recharge vector; {e} T = row vector of l's; and {g} = known right-handside vector reflecting the target gradients. The objective function (3) minimizes the sum of pumping and recharge rates. In this formulation each gradient control well has two separate positive decision variables associated with it, one for pumping and one for recharge. Separating pumping from recharge allows the linear program to pick one or the other, as well as meeting the non-negativity requirement for decision variables. The hydraulic gradient control constraints (4a) guarantee the creation of an inward hydraulic gradient around the shrinking plume perimeter. The target gradient, which is a function o f time, determines the magnitude o f the inward gradient. The formulation of the hydraulic gradient constraints is an extension o f the response matrix method. With this m e t h o d a unit stress is applied at any well location and the drawdown response is determined at any observation location of interest. For each well the groundwater flow model (Trescott et al., 1976) is run for the entire time frame of the problem with a unit stress imposed solely at that well during the first pumping period. At the end o f each pumping period, the drawdown responses are recorded at each observation location. For a meaningful solution of the management problem, the responses generated in the groundwater flow model must be linear to use the principle o f superposition. Superposition allows responses to be multiplied and added together to give the total drawdown response at a given observation point due to all system stresses. As is c o m m o n l y assumed, drawdown is considered

96 to be small compared to the saturated thickness for the unconfined aquifer at the Rocky Mountain Arsenal. The drawdowns are generally within 10-20% o f the saturated thickness in this study. This is necessary to maintain linearity so that superposition and the response matrix technique can be used.

5.3. Derivation o f hydraulic gradient control constraint For this problem the observation locations are the gradient check locations placed around the perimeter of the contaminant plume. The drawdown responses at each pair of two nodes on either side of the plume boundary are converted into a gradient response. In this derivation point A is always considered to be outside the plume and point B is just inside the boundary. The difference in hydraulic head (H) between points A and B should be equal to or greater than a target head difference ( A H t a r g e t ) between the same points. The specified target gradient depends on the problem to be solved; zero for a fiat gradient or a set value to create the desired inward slope: H a - - H b /> A H target

For a system initially at steady state, the head at a given point can be expressed as the steady-state head (/-/steady) minus the drawdown (S) due to stresses in the system at t h a t time: (gsteady.a - - S a ) - - (Hsteady.b - - S b ) ~

Agtarget

Rearranging with known terms on the right gives:

- - ( S ~ - Sb) ~ > - (Hsteady-a --Hsteady-b) + AHtazget We multiply by --1 and define:

Hsteacly.a- -

Hsteady.b

=

AHsteady

This gives the relation for a single pair of points: Sa -- Sb •

AHsteady --

AHtarget

The total drawdown (S) at a given point can be expressed as the sum of the c o m p o n e n t drawdowns due to each of the active pumping and recharge wells. Each c o m p o n e n t drawdown can be expressed as the drawdown due to a unit pumpage at t h a t well, multiplied by the actual well rate. Likewise, the total drawdown (S) due to all pumping and recharge wells can be represented using a unit drawdown response vector {ra} multiplied by the pumping/recharge rates at each well in the system {q}. Drawdown at a single point A (Sa) is: Sa = {ra}W{q} where {r, }W is a row vector of unit drawdown responses at point A at a given time due to each potential well and {q} is a vector of pumping/recharge rates

97 for each potential well. This can be expanded for the set of all points A outside the plume boundary over all time: {Sa} :

[Ra]{q}

Now {sa } is a vector containing the drawdowns for all the exterior points over time. [Ra] is the matrix of unit drawdown responses due to all the potential wells for all the exterior points over all times. Expressing drawdown at the inner gradient check points, B, in a similar manner and defining the right-hand-side gradient target vector, {g}, gives the following expression for the entire system of potential wells and all the gradient check points: IRa] {q}

-

-

[Rb] {q} ~< {g}

or equivalently: ([Ra] -- [Rb]){q} ~< {g} Define [ R a ] - - [ R b ] = [R] where [R] is a single "gradient" response matrix with each coefficient being the differences in unit drawdown responses between adjacent points A and B (gradient check pair): [R]{q} ~< {g} Substitute {q} = { u } - {v}, where {u} = non-negative pumping vector and {v} = non-negative recharge vector. This allows the linear program to pick either pumping {u} or recharge {v} for each well. Minimization of the objective function forces one or the other to always be zero. The final form o f the hydraulic gradient constraint is: [ R ] { u } - - [R]{v}~<{g} The management model has large time divisions called management periods with equal length subdivisions called pumping periods. Each management period is associated with a given position o f the plume during which the number o f potentially active well sites and the hydraulic gradient check locations do n o t change. The pumping/recharge rates at each of the wells can change from pumping period to pumping period, b u t no new potential well sites can be added until the following management period. The reason for the pumping periods is to allow the stress rates to change within the longer management periods. In this problem, the length of the pumping period is the maximum amount o f time allowed to pass before meeting the hydraulic gradient constraints. This forces the constraints to be initially met early in the management period. 6. FORMULATION FOR THE ROCKY MOUNTAIN ARSENAL EXAMPLE Table I shows the specific design o f the R o c k y Mountain Arsenal management problem. The management horizon of 16 yr. is divided into two 8-yr.

98 TABLE I Rocky Mountain Arsenal management problem Management period I

Management period II

Total

Design components: Planning horizon (yr.) Pumping periods Gradient control pairs Potentially active wells

0--8 16 18 9

8--16 16 14 15

16 32 32 15

Matrix statistics: Rows Columns Matrix elements

288 320

224 512

512 832 198,258

management periods, each associated with a plume position that was determined in stage I {Fig. 8). The management periods are subdivided into sixteen ~-yr. pumping periods for a total of 32 pumping periods. There are nine potentially active gradient control wells in management period I. Six more well sites are added in management period II for a total of fifteen active sites. The central contaminant removal well pumps at a fixed rate t h r o u g h o u t the problem. It is assumed that there is an existing network of wells at the beginning of the problem. This means the cost of establishing a well is not taken into account in the optimization problem. It does not mean, however, that the locations of wells are arbitrary and do not affect the solution of the problem. The locations of the wells are very important for obtaining a feasible solution. During trial runs of the linear program, poor well site selection produced infeasible solutions. All the hydraulic gradient constraints could not be met. The infeasibilities were overcome with better selection of well sites or adding new well sites based on the location of the infeasibilities. While not shown here the problem of only selecting a certain number of wells can be handled by using mixed-integer programming. The target gradient specified in the hydraulic constraint determines the magnitude of the inward gradient at the perimeter of the plume. For this problem, the gradient must be at least as large as the gradient produced by pumping the contaminant removal well in an initially fiat hydraulic surface, i.e. the stage-I simulation. With time, the drawdown due to the central well becomes greater and the cone of depression grows. This means the target gradient at the perimeter o f the plume becomes greater with time. The groundwater flow model (Trescott et al., 1976) was run for each of the fifteen gradient control wells with a unit stress of 2.83 1 s -1. The drawdown responses generated in these simulations were converted into gradient

99

mlllll

m

m

m

i

iiiiimll

Iiiiiiii

i

mlllllll

m

'-~ ,/ ~J / /t

~ 0

lnitial Plume Boundary ~ 0

~.,, -'

,

,,"

\

", / 4 - °

O

"' ;" .%-;', ~2

N

~~----~ 0

I

Plume ~?BOUndaW aqi_'!i:; ,',: 8 Years

0

500 Meters ~

I I I I I

~

':('""";; No-Flow Boundary • Central Contaminant Removal Well O Well Locations Gradient Check Location Fig. 8. R o c k y M o u n t a i n Arsenal m a n a g e m e n t p r o b l e m design.

responses at each gradient check location for each pumping period using the m e t h o d described above. The resulting constraint matrix contains 832 decision variables (columns) and 512 constraints (rows). We found during the development of this m e t h o d t h a t the convergence criterion for the numerical scheme used in the groundwater flow model was important when the constraint matrix reaches these dimensions. When adding hundreds of numbers together the convergence criterion of 0.001 typically used in finite-difference models does not produce accurate results. For this study a convergence criterion o f 1 • 10-9 was used. The linear programming problem (3) with constraints (4a) and (4b) was solved using the optimization package, MINOS (Murtagh and Saunders, 1980). Computation times between 43.45 and 217.01 min. on a Prime ® 750 computer were required, depending on the

100

starting point for the solution of the linear program. This is roughly equivalent to between 3 and 14.5 min. on an I B M ® 3 0 3 3 mainframe computer.

7. R E S U L T S

The best well selection and their respective pumping and recharge schedules for containing the contaminant plume during aquifer cleanup were determined for the R o c k y Mountain Arsenal example. The two-stage procedure allows a single global optimization to determine all pumping and recharge schedules for all 32 pumping periods. This capability is possible because plume simulation during stage I predicts plume geometry at all future times. Therefore, the potential well sites (over time and space) and gradient check locations can be fixed at the beginning of the problem. The problem is solved for the entire planning horizon in just one run of the simulation--management model. Another solution strategy is to use a series of sequential optimizations for each pumping period. In this strategy, total pumping and recharge rates are minimized for the first 6-month pumping period. The head distribution resulting from this solution is used as the initial condition for optimization during the next 6-month pumping period. This sequential process is repeated for each of the 32 pumping periods. Stage I, which was an essential component of the global solution, is not necessary in the sequential strategy since each optimization is based only upon the past and does not need any future information. The two above solution strategies, global and sequential, are compared in Fig. 9. A sequential strategy would entail periodically updating the plume configuration before determining new well sites. Although this purely sequential formulation might have resulted in different potential well locations for management period II, identical potential wells are used for an accurate comparison. The results from the global optimization show several interesting results. In management period I all the permitted wells are used to control the gradient. The t w o wells on the upgradient end of the plume pump water o u t o f the aquifer. Seven downgradient wells recharge water. The recharge rates are relatively high during the first ~-yr. pumping period to initially reverse the hydraulic gradient. Then the rates drop during the following 1 to 2 yr., indicating that somewhat lower recharge rates are suitable to maintain the inward gradient. During the remaining portion o f management period I there is a slow rise in the recharge rates. This rise is due to the increase in the target gradient with time. Recall that the target gradient at the plume perimeter becomes greater as the cone of depression around the center well grows. Increased recharge rates are necessary to meet the increasing target gradient. The pumping wells show a similar trend to the recharge wells. Initially, the rates are quite high and in 2--3 yr. the rates level off.

101

SEQUENTIAL SOLUTION

GLOBAL SOLUTION

~

IHHHLUJ

Period 1

Period [

v: 1

"

~: E3 z < 7~

2

v = , i 0 2

i

, = , 4 6 YEAR

i

i 8

--

2"53 I ~

~o.oj

-2.5

2o~

Pumping 1

.......

[]

Recharge ~ ~ Yrs.

~

o

~ y~O

~

• Central Contaminant Removal Well • Active Gradient Control Wells Period [! O Inactive Well Locations ~- 2 ! ~ C o n t a m i n a n t Plume J~ 1

Period II

lo

;

o

~-

8

j

10

12 14 YEAR

16

o~-

8

10

12 14 YEAR

16

/ °

Fig. 9. Global and sequential optimal pumping/recharge schedules for active hydraulic gradient control wells with contaminant removal well pumping at constant rate (2.83 l s-Z).

102

The results for management period II in the global optimization differ from the first period. O f the fifteen possible well sites, the linear program only chooses nine. Of those nine, six are the new well sites that are closest to the 8-yr. plume boundary. The pumping wells initially have low rates that increase b y "~ 10--20% over the next few pumping periods to steady values. The trends can n o t be generalized for the recharge wells. In the global schedules for management period II the pumping and recharge rates dramatically increase or decrease in the very last pumping period. Given the complexity o f the system we remain uncertain as to the explanation for this behavior. In all o f our testing, the same end-period effect occurred when the optimization was run for a shorter period. A terminal value function could be placed in the objective function to eliminate this behavior (Cummings and McFarland, 1974). Such a function would somehow take into account the impacts on the system occurring after the end of the planning horizon. For this case, we suspect that the pumping wells are shut off because their primary benefit for hydraulic gradient control is in future pumping periods. Pumping wells are a greater distance from the critical area for hydraulic gradient control on the north side of the plume. Therefore, their hydraulic influence lags in time behind that of the recharge wells. Under these circumstances, the most efficient strategy is to use the recharge wells alone to control the hydraulic gradient. A similar pattern in which only recharge wells are used is seen in the sequential solution (discussed later) where future pumping periods are not considered by the model. Without the pumping wells, the recharge well rates readjust. Even with the large readjustments, the total pumping/recharge rate for the last pumping period is ~ 5% less than for the previous periods. Projecting the rates from the previous periods into the last period results in a higher (worse) objective value. The sequential management strategy shows some significant differences from the global optimization. Although the total pumping/recharge rates are comparable for the two optimization strategies, the well selection and pumping schedules are very different. The most striking difference is that pumping is all b u t eliminated. During management period I the hydraulic gradient is completely controlled by the recharge wells alone. Only one pumping well is used during management period II at one-third the rate of the same well in the global opimization. The cumulative pumping/recharge rates (Fig. 10), show that the global solution is ~ 10% better than the sequential solution over the entire planning horizon. During management period I the values are nearly the same, the sequential solution being slightly better (lower). The discrepancy increases between the global and sequential cases during management period II. The four graphs near the center of Fig. 9 show the details of these trends. In the first half o f managment period I, the total pumping and recharge rates are higher in the global solution. By management period II, however, not

103

700

1

1

I

600 iii

Sequential Solution

< ~c 500 uJ C9 rt-

<

"1iJ,,J

Global Solution

400

z

~ 300 w

<

200

c~ 100

0

4

L

I

8 YEARS

12

1t

Fig. 10. Cumulative pumping/recharge rates for the global and sequential optimization strategies.

only are the global rates much lower than the sequential rates b u t they are declining while the sequential rates continue to increase. The above results could be expected. The global optimization strategy has the advantage o f solving the problem over the entire time horizon, and therefore is able to take into account constraints at all times. Foresight is built into the solution. The sequential optimization strategy can only consider the constraints that are in effect during the current pumping period. Although global optimization produces a different and better solution than the sequential strategy, the total rates are quite similar and arguments could be made for using either scheme based on other criteria, such as economic or social considerations.

104

8. MANAGEMENT MODEL VERIFICATION

The results from the global strategy were checked by running the solute transport model (Konikow and Bredehoeft, 1978} with the optimal pumping/recharge rates. Successful hydraulic gradient control will prevent the plume from migrating beyond its original boundaries. Redistribution of the contaminant using optimal hydraulic gradient control is compared to redistribution with no control, as shown in Fig. 11. When there are no hydraulic gradient control wells, pumping the contaminant removal well at the specified constant rate has little effect on aquifer restoration or the migration of the plume. In 4 yr. the plume has started to move outside the model area. After 12yr. the plume has escaped from the system. In contrast, with hydraulic gradient control the contaminant plume is contained within its original boundaries throughout the 15yr. necessary for aquifer restoration. The plume geometry resulting from the optimal rates was also compared to the "ideal" plume geometry (Fig. 7) generated in the stage-I simulation. In general, the plume boundary contours in the global optimal solution are much farther north than in the "ideal" case, indicating that the contaminant is being removed less effectively from the northern portion of the system. This is an apparent consequence of imperfect screening of the regional

®

®

Ij I

I

00

6

0 --

T 0

2

8~

o

:'/"

o i

i.-;._~

I._L_.

0 I

500 Meters I.--. I

I

I

I

I

~

I.--.

Central Contaminant Removal Well

' i ' " ,'-' N o - F l o w Boundary I ~ Constant-Head Boundary Flux Boundary

~\,

0

500 Meters I.--,

I I i i I I

Eoio

~(i

0 t4

f

Well Locations Plume Boundary Contours b y Year

Fig. 11. Plume boundary through time: (A) with hydraulic gradient control; and (B) with no hydraulic gradient control.

105

hydraulic gradient. Even though the plume geometries are different, the predicted plume pattern from stage I is suitable for selecting potential well sites. The criterion for placing the potential hydraulic gradient control wells was that they must be outside the predicted plume boundary so that only the central well would remove contaminated water. Although two recharge wells are within the actual plume at the beginning of the second management period, these wells do not remove contaminated water. In fact these wells tend to dilute contaminated water with clean water at the plume margin. The design criteria are not violated.

9. D I S C U S S I O N A N D C O N C L U S I O N S

This study has developed a procedure for aquifer restoration which designs system specific hydraulic gradient control schemes. The R o c k y Mountain Arsenal example demonstrates this procedure for a field-scale transient problem with heterogeneous aquifer properties and complex geometry. The procedure is general and can be applied to contaminant removal in other situations. Although we chose to minimize total pumping, a linear objective -- pumping costs -- could be considered with a quadratic objective function. A mixed-integer formulation would allow the choice of a fixed number of wells. Preliminary investigations have shown that it is possible to include mass-balance constraints which require the total pumping rate to equal the total recharge rate. Many different contaminant removal systems could be used as long as the positions and rates of the clean-up wells are known over the time frame of the management problem. Alternatively the clean-up rate could be treated parametrically in our procedure. Finally, other contaminant geometries can be used with this procedure. For example, the position o f a contaminant front might be o f interest in a salt-water intrusion control problem. This work contributes to m e t h o d o l o g y in two ways: the manner in which the problem is linearized and the extension of the concept o f drawdown responses to hydraulic gradient responses. The nonlinearities in the system are sidestepped b y dividing the problem into two parts. Contaminant transport simulation is separated from the optimization portion of the problem so that a solution can be achieved using optimization techniques for linear systems. Contaminant transport is simulated first so that well site options are based on the expected geometry o f the contaminant plume. This two-stage approach avoids the nonlinearities that would otherwise prohibit the use of a linear optimization technique. The solution for the entire planning horizon, achieved with global optimization differs significantly from the sequential solution in well selection and individual well rates, although the total accumulated pumping plus recharge rates are quite similar. Other economic or social criteria could make either option more desirable. F o r example, if the cost of imported water for

106

recharge is high, the global solution would be more economical. On the other hand, if initial well costs are the deciding factor, the sequential solution might be the better choice. ACKNOWLEDGEMENTS

We gratefully acknowledge the fine review comments provided b y Irwin R e m s o n and Eric Reichard. The use of computer and software brand names in this report is for identification purposes only and does not imply endorsement b y the U.S. Geological Survey. REFERENCES Aguado, E. and Remson, I., 1974. Groundwater hydraulics in aquifer management. J. Hydraul. Div., Proc. Am. Soc. Cir. Eng., 100(HY1): 103--118. Cumming, R.G. and McFarland, J.W., 1974. Groundwater management and salinity control. Water Resour. Res., 10(5): 909--915. Deininger, R.A., 1970. Systems analysis of water supply systems. Water Resour. Bull., 6(4): 573--579. Gorelick, S.M., 1983. A review of distributed parameter groundwater management modeling methods. Water Resour. Res., 19(2): 305--319. Konikow, L.F., 1975. Hydrogeological maps o f the alluvial aquifer in and adjacent to the Rocky Mountain Arsenal, Colorado. U.S. Geol. Surv., Open-File Rep. 74-342. Konikow, L.F., 1977. Modeling chloride movement in the alluvial aquifer at the Rocky Mountain Arsenal, Colorado. U.S. Geol. Surv., Water-Supply Pap. 2044, 43 pp. Konikow, L.F. and Bredehoeft, J.D., 1978. Computer model of two-dimensional solute transport and dispersion in ground water. In: Techniques of Water-Resources Investit a t i o n s o f the USGS, Book 7, Ch. C2. U.S. Government Printing Office, Washington, D.C., 90 pp. Lee, A.S. and Aronofsky, J.S., 1958. A linear programming model for scheduling crude oil production. J. Pet. Technol., 10(7): 51--54. Maddock III, T., 1972. Algebraic technological function from a simulation model. Water Resour. Res., 8(1): 129--134. Molz, F.J. and Bell, L.C., 1977. Head gradient control in aquifers used for fluid storage. Water Resour. Res., 13(6): 795--798. Murtagh, B.A. and Saunders, M.A., 1980. MINOS/Augmented user's manual. Stanford Univ. Syst. Optim. Lab, Stanford, Calif., Tech. Rep. 80-14, 51 pp. R e m s o n , I. and Gorelick, S.M., 1980. Management models incorporating groundwater variables. In: D. Yaron and C.S. Tapiero (Editors), Operations Research in Agriculture and Water Resources. North-Holland Publishing CO., Amsterdam, pp. 333--356. Rosenwald, G.W. and Green, D.W., 1974. A method for determining the optimum location o f wells in a reservoir using mixed-integer programming. Soc. Pet. Eng. J., 14: 44--54. Scheidegger, A.E., 1961. General theory of dispersion in porous media. J. Geophys. Res., 66(10): 3273--3278. Trescott, P.C., Pinder, G.F. and Larson, S.P., 1976. Finite-difference model for aquifer simulation in two dimensions with results of numerical experiments. In: Techniques of Water-Resources Investigations of the USGS, Book 7, Ch. C1. U.S. Government Printing Office, Washington, D.C., 116 pp. Warner, J.W., 1979. Digital-transport model study of diisopropyl-methyphosphonate (DIMP) groundwater contamination at the R o c k y Mountain Arsenal, Colorado. U.S. Geol. Surv., Open-File Rep. 79-676, 39 pp. Wattenbarger, R.A., 1970. Maximizing seasonal withdrawals from gas storage reservoirs. J. Pet. Technol., 22(8): 994--998.