- Email: [email protected]
Analysis of the impact of layered soil heterogeneity on optimal policies for groundwater, remediation
1159
A n a l y s i s of t h e i m p a c t of l a y e r e d soil h e t e r o g e n e i t y o n o p t i m a l p o l i c i e s for groundwater
remediation
D. Bati a and A.S. Mayer a* aDepartment of Geological and Mining Engineering and Sciences, Michigan Technological University, 1400 Townsend Drive, Houghton MI 49931, USA Typical groundwater remediation problems involve selection of the number, location, and flow rate schedule of pumping and injection wells. Simulation models combined with optimization models are used to find optimal designs while considering management objectives and constraints. Due to the high computational effort involved and lack of data, simulation models often are based on simplified 2-D homogeneous hydrogeological settings. The purpose of this work is to investigate how simplifying hypotheses may affect the final optimal relnediation policy. In particular, the analysis addresses the case of a heterogeneous layered aquifer versus a homogeneous one with an equivalent (lumped) hydraulic conductivity. The analysis is carried out by determining the optimal solution for the remediation of the plume in the case of a homogeneous aquifer, and then simulating this optimal policy with the actual layered aquifer. 1. I N T R O D U C T I O N
The remediation of dissolved chemicals from shallow aquifers is often tackled with pump-and-treat (PAT) techniques. Typical PAT problems involve selection of the number, location, and flow rate schedule of pumping and injection wells, and selection of the most appropriate treatment method. For a given treatment technique, contaminant flow and transport models combined with optimization algorithms are used to find optireal designs while considering management objectives and constraints. Due to high costs for aquifer characterization and to the heavy computational burden involved, simulation models often are based on simplified 2-D homogeneous hydrogeological settings. Hydrogeological measurements, such as hydraulic conductivity, obtained from in situ tests are usually unable to describe the heterogeneity of the porous medium. As long as expensive aquifer characterization is not undertaken, assumptions of homogeneity appear to be the most reasonable choice. The purpose of this work is to investigate how simplifying hypotheses may affect the final optimal remediation policy. The analysis compares a heterogeneous unconfined lay*The authors wish to thank Mario Putti, from the University of Padua, Italy, for the gracious use of his numerical groundwater flow simulator. The majority of this work has been made possible by an EPA Office and Research Development grant (CR826614-01-0) and an NSF grant from the Bioengineeringand Environmental Systems Division (BES-0083112).
1160 ered aquifer to a homogeneous one with an equivalent (lumped) hydraulic conductivity. The example being investigated refers to the benchmark problems outlined by [8] for use in the research community as a means for comparing optimization approaches. Groundwater flow and contaminant transport are simulated with a fully 3-D finite element (FE) unsaturated flow code along with a random walk (RW) particle tracking transport code. The flow and transport code is coupled to a genetic algorithm (GA) [10] code to solve the optimization problem represented by the minimization of the cost of remediating a hypothetical aquifer-contaminant system using PAT. The ez situ treatment technology for the contaminated groundwater is air stripping. The cost function is a nonlinear function of decision variables (pumping rates) and state variables (hydraulic heads and contaminant concentrations). Constraints include limits on hydraulic head and the contaminant mass remaining in the aquifer at the end of the remedial horizon. The results of homogeneous and heterogeneous simulations are compared in terms of cost and values of the decision variables: the hypothetical plume is optimally cleaned up as though only an equivalent value of hydraulic conductivity was known, i.e., measured by means of in situ pumping tests. The optimal policy, in the form of a set of pumping rates in as many pumping wells, is then applied to the actual layered aquifer. The two solutions are compared in terms of cleanup performance. 2.
MODEL
2.1.
EQUATIONS
Variably
saturated
flow
Isothermal fluid flow in variably saturated 3-D porous media is governed by the partial differential equation (PDE) obtained by substituting Darcy's law into the continuity equation [1]. Using the pressure head ~ as the dependent variable the PDE reads: d~/
--~
q ; i,j-1,2,3
(1)
where: -
xi
t
-
reference system coordinates (L); time(T); saturated hydraulic conductivity tensor (L/T);
=
relative conductivity with respect to water phase (dimensionless);
-
saturation of water (dimensionless);
S~ q f]l
specific elastic storage (l/L);
-
porosity (dimensionless);
-
source/sink (per unit volume) (I/T);
--
f/2 --
1;f/3
--
0.
Equation (1) is nonlinear due to pressure head dependencies in the storage and conductivity terms. In this work no dependency of S~ on ~ is considered, whereas saturation S~ and relative conductivity K ~ are expressed as functions of pressure head using the characteristic relations from [2].
1161 2.2. F i n i t e e l e m e n t s o l u t i o n
PDE (1) is solved using a Galerkin FE discretization in space and linear basis functions, with 3-D tetrahedrons. A weighted finite difference scheme is adopted for time discretization [3]. The ensuing nonlinear set of equations may be linearized and solved using the following Picard iteration [4]"
[
pH(, j 1 plk~+))]. (k+,) + At(k)
~)(k+L,) -
E1
(k+l)
P ( k + , ) - (1
-('~)
At(k)
/J~/)(k+l) +
1 (1-- ") ~(k) ,9 -2- <
-<
i
(2)
where (k) is the time-step index, (m) is the iteration index, and H and P are the flow stiffness and capacity matrices, respectively. The velocity field may be determined using Darcy's law from the solution of (2). The component of pore velocity along xi is then given by:
+ where 0 - r
(3)
is the soil moisture content.
2.3. T r a n s p o r t of n o n - r e a c t i v e
solute
The transport of conservative solutes in a partially saturated porous medium relies on the PDE [5]"
0 [ O~xj] O(viC) OC Oxi Dij - Oxi = Ot
qC*+f, 0
'
i j-123 ' ' '
(4)
where C is the solute concentration (M/L3), C* is the concentration associated with q (M/L3), and f represents the specific solute mass released with no fluid exchange (ML3/T). Dij, the classic hydrodynamic dispersion tensor (L2/T), is written as [5]" ViVj
D*
(5)
where 6id is the Kronecker delta function, Iv] - I E v ~ (L/T), C~L and aT are the longitudinal and transverse dispersivities (L), T is the tortuosity; and D* is the free liquid diffusivity (L2/T). 2.4. P a r t i c l e t r a c k i n g s o l u t i o n
The transport PDE (4) is addressed with a particle tracking method based on a RW algorithm [6]. Following this technique, the solute is represented as a large collection of particles each of which is assigned a mass. For a given time step, each particle is displaced in two separate stages. The first involves a translation along a streamline computed from the pore velocity distribution (Equation (3)); the second involves a random displacement, the direction and magnitude of which are determined so that the overall distribution of particles will resemble the solution of (4). The RW solution for spatially dependent
1162 dispersive transport under unsaturated conditions reads [7]: X}k-I-l)
--
1 D}k) c90 (k)) /~tk -Jr- blk)Zj ~///~k X}k)-71- V}k)~ OD}~ OXj ) Jr--~k OXj 1
D}~ ) - -~bl~)b~ ~ ; Zj - N(O, 1)
(6) (7)
where x~k) indicates coordinate xi of the generic particle at time tk. The third term at the right-hand side of (6) represents the random nature of hydrodynamic dispersion. In Equation (7) vector Zj (j - 1,2,3) is obtained from three independent random numbers picked from a standard normal distribution N(0, 1). In this study, we elected to neglect the two terms at the right-hand side of (6) that account for the spatial gradients of hydrodynamic dispersion and soil moisture content. 3. D E S C R I P T I O N
OF GROUNDWATER
SYSTEM
3.1. O v e r v i e w
In the simulated test case the domain is represented by a 10001000 (ram), 30-m thick unconfined aquifer, wherein regional steady-state flow is determined by the boundary conditions represented in Figure 1. The aquifer extents between z=0 in and z=30 m and all elevations are taken with respect to the aquifer bottom (z=0 m). The bottom of the aquifer is impermeable, and the aquifer is recharged with a uniform and constant 1.28210 -5 m/s rate. The solute plume develops from a 25502 source (Figure 1) just below the groundwater table within which a concentration C~-1 kg/m 3 is enforced to persist for a period of 5 years and removed thereafter. The aquifer heterogeneity is represented by 10 equal-thickness homogeneous layers with variable hydraulic conductivity as shown in Figure 2. The equivalent horizontal hydraulic conductivity is given by the arithmetic mean of these 10 values and is equal to 5.0110 -5 m/s. Other hydrogeological parameters can be found in Table 1; more details can be found in [8]. 3.2. N u m e r i c a l t e s t s
Numerical tests have been performed with the hydrogeological parameters listed in Table 1. Figure 3a shows the horizontal projection of the 3-D FE grid used to simulate the regional steady-state flow field over which the contaminant is released. The corresponding water table position is contoured in Figure 4 for both the homogeneous and the layered
Table 1 Hydrogeological parameters for the simulated test case. Parameter Symbol Value Units 5.0110 -5 m/s Average Hydraulic Conductivity K0 Porosity r 0.32 1 m Longitudinal Dispersivity ctc 0.1 m Transverse Dispersivity aT 10 -9 m2/s/ Molecular Diffusivity D* Tortuosity ~0.4
1163
Figure 1. Flow and transport boundary conditions for the simulated test case.
Figure 2. Distribution of hydraulic conductivity in the 10-layer aquifer. The equivalent horizontal hydraulic conductivity is given by the arithmetic mean of these 10 values and is equal to 5.0110 -~ m/s.
aquifer. It may be observed that the two solutions are quite similar to each other, while more significant discrepancies are to be expected in the local pore velocity due to the high variability of hydraulic conductivity among layers (Figure 2). The pore velocity and the water saturation fields in the layered aquifer as obtained from the FE model (Equations (2-3)) are used in the RW particle tracking code (Equations (47)) to determine the position of the contaminant plume to be cleaned up. Contaminant particles originate from a 25502 source (Figures 1 and 4) just below the groundwater table in order to enforce a uniform and constant concentration C~= 1 kg/m 3 for a period of 5 years in the release region. In the RW simulation a time step A t = l day was applied with 1 particle per kg solute mass. The generated plume, represented in Figure 5 in terms of particle positions, consists of 3.5104 kg of solute mass. It is worth noting that the contaminant tends to concentrate in the low permeable layers above and below the source. Figure 5 also shows the locations chosen for 8 potential pumping wells. The wells are screened for the entire thickness of the aquifer. This condition is modeled by imposing a very large value for the hydraulic conductivity in the elements around the well nodes, so as to maintain a constant hydraulic head throughout the elements in the well itself. Figures 3b, 3c, and 3d display the horizontal projection of the 3-D FE grid assembled to simulate the flow field when pumping wells are activated. Wells are set up along the center of the plume approximately every 50 meters. Each contaminant particle is removed from the system when it is sufficiently close to the pumping well position, in this case at a distance of about 2 meters (Figure 3d). This condition is imposed in order
1164
Figure 3. Horizontal projections of the FE grid used to simulate" (a) regional steady-state flow; (b) flow field with 8 potential pumping wells activated for remediation. Particulars of mesh (b)are also shown ((c)and (d)).
to overcome particle overshooting effects due to inaccurate pore velocities determined in elements connected to source/sink regions. In Figure 6, the average concentration in some of the layers below the groundwater table is plotted. As already mentioned (Figure 5), it may be observed that, due to hydrodynamic dispersion, most of the contaminant ended up in the low permeability layers below the water table (Figure 2), and, to a lesser extent, above in the vadose zone. As matter of fact, we are aware that the concentration calculated in these "stagnant" regions may suffer from inaccuracies since dispersion and soil moisture gradients in Equation (6) have been overlooked [9]. However, the cleanup of low permeable layers of the aquifer is expected to require both high pumping rates and a long operation time. By consequence, the optimal PAT policy is expected to be strongly influenced by the assumption of either a more realistic layered aquifer or an equivalent homogeneous one. In what follows the plume "created" in the layered aquifer will be optimally cleaned up as though an equivalent value of hydraulic conductivity was known, for example, as would be measured by means of in situ pumping tests. The optimal set of pumping rates in
1165
Figure 4. Contour representation of the water table position in the regional steady-state flow assuming both a homogeneous (solid lines) and a layered (dashed lines) aquifer.
Figure 5. Particle positions in the plume developed after 5 years of contaminant release. The locations of 8 potential wells are also shown.
the potential 8 pumping wells is then applied to the "real" layered aquifer. No injection wells will be considered in this case. The two solutions are compared in terms of cleanup performance. 4. M A N A G E M E N T
PROBLEM
A single-objective optimization algorithm is applied to the remediation of the contaminated aquifer by PAT. The optimization framework uses the Pareto genetic algorithm (GA) [10] and is applied to minimize the objective function while requiring a series of constraints to be met. The decision variables are extraction rates at 8 potential pumping well locations, shown in Figure 5. During the considered operation period pumping rates are constant and steady-state flow is assumed to be reached instantaneously at the start of remediation. The ez situ treatment technology for the contaminated groundwater is air stripping. 4.1. O b j e c t i v e function The objective function (OF) to be minimized is represented by the cost of the PAT remediation system, which includes capital costs incurred from installation of recovery wells and operational costs from pumping and groundwater treatment [8]"
OF Total Cost Capital Cost
9 -
min{Total
Capital Cost + Operational
:
(8)
Cost} Cost - F(Q~)
i -
1, 2, .., n
Well Construction, Pumps, and Stripping Tower costs
(9)
1166
Figure 6. Average concentration in some of the layers below the groundwater table at the start of remediation (t - 5 years). The majority of the solute mass is located in the low permeable layer between the elevations of 15 and 18 m. ?%
Capital Cost
Tte
E cod~~ --~ E c1 [QT[ bl (zg~- h.~i~) b2 + c4Z
-
i=1
Operational Cost
9 Well and P u m p Operation, and Treatment costs
Operational Cost
-
Jtin
(10)
i=1
c2Q~(h~ - zg8) -}- E c3Qi + Z(c5 - C6QT) dt i=ne +1
(11) where" n
=
total number of injection/extraction wells;
~e
=
number of extraction wells;
Qi
injection/extraction rate for well i (m3/s);
-
-
hmin Zgs
di tin, t f in
Z
Z
design extraction rate for well i (m3/s);
-
total extraction rate (m 3/s);
-
hydraulic head in well i (m); minimum allowable hydraulic head (m);
-
-
-
-
elevation of ground surface (m); depth of well i below the ground surface (m); initial and final time for remediation (s); air-stripping tower height (m);
]Q~[b5
[Ce--CElb6c?}
-
max{b3H b4
=
Henry's law coefficient = 0.2;
,
_> c E ) ;
influent concentration (kg/m3); CE
=
target effluent concentration = 5 10 .5 k g / m 3.
The coefficients bi and ci associated with the OF may be found in [8].
1167 4.2. C o n s t r a i n t s The state and decision variables are constrained by both hydraulic and solute mass constraints [8]. Hydraulic constraints include limitations on the extraction and injection rates at each well"
o ~e xa x <- -Q i < o- -m a xi n
i -- 1, 2, .., n
(12)
n
OT~ X
< - - ~-~Oi
(13)
i=1
In (12) ~ 0 (~max
max <
and -~i~ 0 "~x are the maximum extraction and injection rates at any well
0
--
" (~max ~
>
_
0
" [(~max
I __ Io m a x I - 6.410 "* in
-3
I113
/s
In (13) I O m a x is the maximum total extraction rate which has been set to -3 210 .2 m3/s Constraints on pumping rates (Equations (12)-(13)) are established by setting lower and upper bounds to the values of the decision variables between which the GA searches for the optimum. Hydraulic constraints include also limitations on the aquifer drawdown" '~T
~
hmin <__hi < hmax i -
"
"
(14)
1, 2, .., n
where h,~i~- 10 m and hmax=Zgs. Solute mass constraints are specified by requiring (a) the contaminant concentration to be zero on the fixed head boundaries shown in Figure 1 at any time during remediation; (b) a minimum aquifer cleanup performance as quantified by the global fraction of mass that may remain in the system at the end of the remediation horizon" C(x-
1000, y, z, t ~ < t <_ t f ~ ) -
C(x, y -
1000, z, t ~ < t <_ t f ~ ) -
0
(15) (16)
with the operation time t fin-tin equal to 10 years, and the PAT goal M set to 0.05. Both drawdown (Equation (14)) and solute mass constraints (Equations (15)-(16)) are enforced with high-order-of-magnitude penalty coefficients that simply multiply the total cost function if the constraints are not met. 5. R E S U L T S
AND
DISCUSSION
The management problem is solved by linking a GA code [10,11] to the F E / R W flow and transport model described in Section 2, which represents the "function call" calculating the cost associated with each particular design. In the equivalent-K (homogeneous) test case the optimum was found with about 470 OF evaluations. Figure 7 shows the optimal solution determined by the GA. Only 4 of the 8 potential wells are included in the optimal design. The maximum flow rate at well 8 is needed in order to prevent contaminant from reaching the upper boundary of the domain as prescribed by the constraint of Equation (15). Lower flow rates are needed from the other wells in the higher concentration areas of the plume so as to satisfy the minimum cleanup performance as in Equation (16).
1168
Figure 7. GA results for the 8-well design: Q2=Q6=4.63010-4m3/s (40 m3/day); Q4=2.31510-4m3/s (20 m3/day); Qs-3.47210-3m3/s (300 ran/day); wells 1, 3, 5, 7 are inactive.
Figure 8. Distribution of contaminant particles in the layered aquifer at the end of the 10-years long remedial time, where optimal policy obtained from equivalentK (homogeneous) case is applied to the layered case.
Table 2 reports the costs and the solute mass remaining associated with this solution. Table 2 also reports the results obtained when the optimal remediation policy determined with the equivalent-K (homogeneous) aquifer is applied to the heterogeneous aquifer. As anticipated, the performance of the policy calculated for the equivalent homogeneous aquifer is very poor. About 32% of the solute mass remains in the aquifer at the end of the remedial horizon. Figure 8 describes the distribution of contaminant particles in the layered aquifer 10 years after operations began. Solute remains primarily in the low permeable layer between the elevations of 15 and 18 In and the flow rates obtained from the equivalent-K case are not adequate.
Table 2 Comparison of remediation costs and performances. Cost ($) Aquifer Mr (%) Capital Pump Treat Homogeneous 0.5 182,380 6,790 63,770 Layered 32.6 185,630 6,960 65,740
Total 252,940 258,330(*)
(*) The cost does not reflect the p e n a l t y to be applied because constraint (16) is not met.
1169
Figure 9. Total mass of solute (kg) in the aquifer vs time. The period of time spans over 15 years from the start of contaminant release up until the end of the remedial operations.
In Figure 9, the behaviors of the total solute mass in the aquifer vs. time are displayed. The solid profile represents the plume development as simulated for the heterogeneous aquifer. The dashed profile indicates how contaminant is extracted from the aquifer following the optimal solution (Figure 7) determined assuming an equivalent homogeneous porous medium. The substantially double-slope shape of the diagram suggests that a solution with variable pumping rates would probably allow reducing the total cost. Finally, the dotted profile reveals how the contaminant would be cleaned up if the optimal policy of Figure 7 were applied to the layered heterogeneous aquifer. These results, though obtained using two simplistic conceptual models for the aquifer, the homogenous and the layered one, prove how important site characterization and collection of data may be in order to determine the optimal solution for the PAT process. We note, however, that the quantitative results presented here may be somewhat biased by neglecting the dispersion and soil moisture gradient terms in Equation (6).
1170 6. C O N C L U S I O N S
In this work we analysed how simplifying hypotheses may affect the final optimal remediation policy. A heterogeneous unconfined layered aquifer was compared to a homogeneous one with an equivalent (lumped) hydraulic conductivity. An optimal design was obtained for the homogeneous aquifer and this design was applied to the layered aquifer. Analysis of the two solutions, compared in terms of cleanup performance associated with the realization of the process, revealed that the performance of the policy calculated for the equivalent homogeneous aquifer results was in fact very poor. At the end of the planned remedial horizon, more than 30% of the solute is still in the aquifer, located in low permeable layers below the water table. The cleanup of these ~'stagnant" regions of the aquifer thus requires pumping rates and operation times much higher than those estimated in the homogeneous scenario. These results prove the importance of site characterization in order to determine the optimal solution for the PAT process. Future development of this work will consider the inclusion of site characterization in the objective function to evaluate when a more thorough knowledge of the aquifer is justified in term of both cost and remediation efficiency; the presence of both injection and extraction wells; the comparison of results for layered aquifers versus aquifers with higher degrees of heterogeneity; and the inclusion of dispersion and soil moisture gradient terms. REFERENCES
1. J . R . Philip, Adv. Hydrosci., 5 (1969) 215. 2. P.S. Huyakorn, S. D. Thomas, and B. M. Thompson, Water Resour. Res. 20 (1984) 1099. 3. P.S. Huyakorn and G. F. Pinder, Computational Methods in Subsurface Flow, 1983, Academic, San Diego. 4. C. Paniconi and M. Putti, Water Resour. Res. 30 (1994) 3357. 5. J. Bear, Hydraulics of Groundwater, 1979, McGraw-Hill, New York. 6. G. Uffink, NATO Advanced Workshop on Advances in Analytical and Numerical Groundwater Flow and Quality Modeling, Portugal, 1987. 7. A . F . B . Tompson, E. G. Vomvoris and L. W. Gelhar, Numerical Simulation of Solute Transport in Randomly Heterogeneous Porous Media: Motivation, Model development, and Application. Rep. UCID-21281, Lawrence Livermore Natl. Lab., Livermore, California, (1987). 8. A.S. Mayer, C. T. Kelley and C. T. Miller, Adv. Water Resour. 25 (2002) 1233. 9. E . R . LaBolle G. E. Fogg and A. F. B. Tompson, Water Resour. Res. 32 (1996) 583. 10. D. E. Golberg, Genetic algorithms in search, optimization, and machine learning, 1989, Addison-Wesley, Reading, MA. 11. M. Erickson, A. S. Mayer and J. Horn, Adv. Water Resour. 25 (2002) 51.
A n a l y s i s of t h e i m p a c t of l a y e r e d soil h e t e r o g e n e i t y o n o p t i m a l p o l i c i e s for groundwater
remediation
D. Bati a and A.S. Mayer a* aDepartment of Geological and Mining Engineering and Sciences, Michigan Technological University, 1400 Townsend Drive, Houghton MI 49931, USA Typical groundwater remediation problems involve selection of the number, location, and flow rate schedule of pumping and injection wells. Simulation models combined with optimization models are used to find optimal designs while considering management objectives and constraints. Due to the high computational effort involved and lack of data, simulation models often are based on simplified 2-D homogeneous hydrogeological settings. The purpose of this work is to investigate how simplifying hypotheses may affect the final optimal relnediation policy. In particular, the analysis addresses the case of a heterogeneous layered aquifer versus a homogeneous one with an equivalent (lumped) hydraulic conductivity. The analysis is carried out by determining the optimal solution for the remediation of the plume in the case of a homogeneous aquifer, and then simulating this optimal policy with the actual layered aquifer. 1. I N T R O D U C T I O N
The remediation of dissolved chemicals from shallow aquifers is often tackled with pump-and-treat (PAT) techniques. Typical PAT problems involve selection of the number, location, and flow rate schedule of pumping and injection wells, and selection of the most appropriate treatment method. For a given treatment technique, contaminant flow and transport models combined with optimization algorithms are used to find optireal designs while considering management objectives and constraints. Due to high costs for aquifer characterization and to the heavy computational burden involved, simulation models often are based on simplified 2-D homogeneous hydrogeological settings. Hydrogeological measurements, such as hydraulic conductivity, obtained from in situ tests are usually unable to describe the heterogeneity of the porous medium. As long as expensive aquifer characterization is not undertaken, assumptions of homogeneity appear to be the most reasonable choice. The purpose of this work is to investigate how simplifying hypotheses may affect the final optimal remediation policy. The analysis compares a heterogeneous unconfined lay*The authors wish to thank Mario Putti, from the University of Padua, Italy, for the gracious use of his numerical groundwater flow simulator. The majority of this work has been made possible by an EPA Office and Research Development grant (CR826614-01-0) and an NSF grant from the Bioengineeringand Environmental Systems Division (BES-0083112).
1160 ered aquifer to a homogeneous one with an equivalent (lumped) hydraulic conductivity. The example being investigated refers to the benchmark problems outlined by [8] for use in the research community as a means for comparing optimization approaches. Groundwater flow and contaminant transport are simulated with a fully 3-D finite element (FE) unsaturated flow code along with a random walk (RW) particle tracking transport code. The flow and transport code is coupled to a genetic algorithm (GA) [10] code to solve the optimization problem represented by the minimization of the cost of remediating a hypothetical aquifer-contaminant system using PAT. The ez situ treatment technology for the contaminated groundwater is air stripping. The cost function is a nonlinear function of decision variables (pumping rates) and state variables (hydraulic heads and contaminant concentrations). Constraints include limits on hydraulic head and the contaminant mass remaining in the aquifer at the end of the remedial horizon. The results of homogeneous and heterogeneous simulations are compared in terms of cost and values of the decision variables: the hypothetical plume is optimally cleaned up as though only an equivalent value of hydraulic conductivity was known, i.e., measured by means of in situ pumping tests. The optimal policy, in the form of a set of pumping rates in as many pumping wells, is then applied to the actual layered aquifer. The two solutions are compared in terms of cleanup performance. 2.
MODEL
2.1.
EQUATIONS
Variably
saturated
flow
Isothermal fluid flow in variably saturated 3-D porous media is governed by the partial differential equation (PDE) obtained by substituting Darcy's law into the continuity equation [1]. Using the pressure head ~ as the dependent variable the PDE reads: d~/
--~
q ; i,j-1,2,3
(1)
where: -
xi
t
-
reference system coordinates (L); time(T); saturated hydraulic conductivity tensor (L/T);
=
relative conductivity with respect to water phase (dimensionless);
-
saturation of water (dimensionless);
S~ q f]l
specific elastic storage (l/L);
-
porosity (dimensionless);
-
source/sink (per unit volume) (I/T);
--
f/2 --
1;f/3
--
0.
Equation (1) is nonlinear due to pressure head dependencies in the storage and conductivity terms. In this work no dependency of S~ on ~ is considered, whereas saturation S~ and relative conductivity K ~ are expressed as functions of pressure head using the characteristic relations from [2].
1161 2.2. F i n i t e e l e m e n t s o l u t i o n
PDE (1) is solved using a Galerkin FE discretization in space and linear basis functions, with 3-D tetrahedrons. A weighted finite difference scheme is adopted for time discretization [3]. The ensuing nonlinear set of equations may be linearized and solved using the following Picard iteration [4]"
[
pH(, j 1 plk~+))]. (k+,) + At(k)
~)(k+L,) -
E1
(k+l)
P ( k + , ) - (1
-('~)
At(k)
/J~/)(k+l) +
1 (1-- ") ~(k) ,9 -2- <
-<
i
(2)
where (k) is the time-step index, (m) is the iteration index, and H and P are the flow stiffness and capacity matrices, respectively. The velocity field may be determined using Darcy's law from the solution of (2). The component of pore velocity along xi is then given by:
+ where 0 - r
(3)
is the soil moisture content.
2.3. T r a n s p o r t of n o n - r e a c t i v e
solute
The transport of conservative solutes in a partially saturated porous medium relies on the PDE [5]"
0 [ O~xj] O(viC) OC Oxi Dij - Oxi = Ot
qC*+f, 0
'
i j-123 ' ' '
(4)
where C is the solute concentration (M/L3), C* is the concentration associated with q (M/L3), and f represents the specific solute mass released with no fluid exchange (ML3/T). Dij, the classic hydrodynamic dispersion tensor (L2/T), is written as [5]" ViVj
D*
(5)
where 6id is the Kronecker delta function, Iv] - I E v ~ (L/T), C~L and aT are the longitudinal and transverse dispersivities (L), T is the tortuosity; and D* is the free liquid diffusivity (L2/T). 2.4. P a r t i c l e t r a c k i n g s o l u t i o n
The transport PDE (4) is addressed with a particle tracking method based on a RW algorithm [6]. Following this technique, the solute is represented as a large collection of particles each of which is assigned a mass. For a given time step, each particle is displaced in two separate stages. The first involves a translation along a streamline computed from the pore velocity distribution (Equation (3)); the second involves a random displacement, the direction and magnitude of which are determined so that the overall distribution of particles will resemble the solution of (4). The RW solution for spatially dependent
1162 dispersive transport under unsaturated conditions reads [7]: X}k-I-l)
--
1 D}k) c90 (k)) /~tk -Jr- blk)Zj ~///~k X}k)-71- V}k)~ OD}~ OXj ) Jr--~k OXj 1
D}~ ) - -~bl~)b~ ~ ; Zj - N(O, 1)
(6) (7)
where x~k) indicates coordinate xi of the generic particle at time tk. The third term at the right-hand side of (6) represents the random nature of hydrodynamic dispersion. In Equation (7) vector Zj (j - 1,2,3) is obtained from three independent random numbers picked from a standard normal distribution N(0, 1). In this study, we elected to neglect the two terms at the right-hand side of (6) that account for the spatial gradients of hydrodynamic dispersion and soil moisture content. 3. D E S C R I P T I O N
OF GROUNDWATER
SYSTEM
3.1. O v e r v i e w
In the simulated test case the domain is represented by a 10001000 (ram), 30-m thick unconfined aquifer, wherein regional steady-state flow is determined by the boundary conditions represented in Figure 1. The aquifer extents between z=0 in and z=30 m and all elevations are taken with respect to the aquifer bottom (z=0 m). The bottom of the aquifer is impermeable, and the aquifer is recharged with a uniform and constant 1.28210 -5 m/s rate. The solute plume develops from a 25502 source (Figure 1) just below the groundwater table within which a concentration C~-1 kg/m 3 is enforced to persist for a period of 5 years and removed thereafter. The aquifer heterogeneity is represented by 10 equal-thickness homogeneous layers with variable hydraulic conductivity as shown in Figure 2. The equivalent horizontal hydraulic conductivity is given by the arithmetic mean of these 10 values and is equal to 5.0110 -5 m/s. Other hydrogeological parameters can be found in Table 1; more details can be found in [8]. 3.2. N u m e r i c a l t e s t s
Numerical tests have been performed with the hydrogeological parameters listed in Table 1. Figure 3a shows the horizontal projection of the 3-D FE grid used to simulate the regional steady-state flow field over which the contaminant is released. The corresponding water table position is contoured in Figure 4 for both the homogeneous and the layered
Table 1 Hydrogeological parameters for the simulated test case. Parameter Symbol Value Units 5.0110 -5 m/s Average Hydraulic Conductivity K0 Porosity r 0.32 1 m Longitudinal Dispersivity ctc 0.1 m Transverse Dispersivity aT 10 -9 m2/s/ Molecular Diffusivity D* Tortuosity ~0.4
1163
Figure 1. Flow and transport boundary conditions for the simulated test case.
Figure 2. Distribution of hydraulic conductivity in the 10-layer aquifer. The equivalent horizontal hydraulic conductivity is given by the arithmetic mean of these 10 values and is equal to 5.0110 -~ m/s.
aquifer. It may be observed that the two solutions are quite similar to each other, while more significant discrepancies are to be expected in the local pore velocity due to the high variability of hydraulic conductivity among layers (Figure 2). The pore velocity and the water saturation fields in the layered aquifer as obtained from the FE model (Equations (2-3)) are used in the RW particle tracking code (Equations (47)) to determine the position of the contaminant plume to be cleaned up. Contaminant particles originate from a 25502 source (Figures 1 and 4) just below the groundwater table in order to enforce a uniform and constant concentration C~= 1 kg/m 3 for a period of 5 years in the release region. In the RW simulation a time step A t = l day was applied with 1 particle per kg solute mass. The generated plume, represented in Figure 5 in terms of particle positions, consists of 3.5104 kg of solute mass. It is worth noting that the contaminant tends to concentrate in the low permeable layers above and below the source. Figure 5 also shows the locations chosen for 8 potential pumping wells. The wells are screened for the entire thickness of the aquifer. This condition is modeled by imposing a very large value for the hydraulic conductivity in the elements around the well nodes, so as to maintain a constant hydraulic head throughout the elements in the well itself. Figures 3b, 3c, and 3d display the horizontal projection of the 3-D FE grid assembled to simulate the flow field when pumping wells are activated. Wells are set up along the center of the plume approximately every 50 meters. Each contaminant particle is removed from the system when it is sufficiently close to the pumping well position, in this case at a distance of about 2 meters (Figure 3d). This condition is imposed in order
1164
Figure 3. Horizontal projections of the FE grid used to simulate" (a) regional steady-state flow; (b) flow field with 8 potential pumping wells activated for remediation. Particulars of mesh (b)are also shown ((c)and (d)).
to overcome particle overshooting effects due to inaccurate pore velocities determined in elements connected to source/sink regions. In Figure 6, the average concentration in some of the layers below the groundwater table is plotted. As already mentioned (Figure 5), it may be observed that, due to hydrodynamic dispersion, most of the contaminant ended up in the low permeability layers below the water table (Figure 2), and, to a lesser extent, above in the vadose zone. As matter of fact, we are aware that the concentration calculated in these "stagnant" regions may suffer from inaccuracies since dispersion and soil moisture gradients in Equation (6) have been overlooked [9]. However, the cleanup of low permeable layers of the aquifer is expected to require both high pumping rates and a long operation time. By consequence, the optimal PAT policy is expected to be strongly influenced by the assumption of either a more realistic layered aquifer or an equivalent homogeneous one. In what follows the plume "created" in the layered aquifer will be optimally cleaned up as though an equivalent value of hydraulic conductivity was known, for example, as would be measured by means of in situ pumping tests. The optimal set of pumping rates in
1165
Figure 4. Contour representation of the water table position in the regional steady-state flow assuming both a homogeneous (solid lines) and a layered (dashed lines) aquifer.
Figure 5. Particle positions in the plume developed after 5 years of contaminant release. The locations of 8 potential wells are also shown.
the potential 8 pumping wells is then applied to the "real" layered aquifer. No injection wells will be considered in this case. The two solutions are compared in terms of cleanup performance. 4. M A N A G E M E N T
PROBLEM
A single-objective optimization algorithm is applied to the remediation of the contaminated aquifer by PAT. The optimization framework uses the Pareto genetic algorithm (GA) [10] and is applied to minimize the objective function while requiring a series of constraints to be met. The decision variables are extraction rates at 8 potential pumping well locations, shown in Figure 5. During the considered operation period pumping rates are constant and steady-state flow is assumed to be reached instantaneously at the start of remediation. The ez situ treatment technology for the contaminated groundwater is air stripping. 4.1. O b j e c t i v e function The objective function (OF) to be minimized is represented by the cost of the PAT remediation system, which includes capital costs incurred from installation of recovery wells and operational costs from pumping and groundwater treatment [8]"
OF Total Cost Capital Cost
9 -
min{Total
Capital Cost + Operational
:
(8)
Cost} Cost - F(Q~)
i -
1, 2, .., n
Well Construction, Pumps, and Stripping Tower costs
(9)
1166
Figure 6. Average concentration in some of the layers below the groundwater table at the start of remediation (t - 5 years). The majority of the solute mass is located in the low permeable layer between the elevations of 15 and 18 m. ?%
Capital Cost
Tte
E cod~~ --~ E c1 [QT[ bl (zg~- h.~i~) b2 + c4Z
-
i=1
Operational Cost
9 Well and P u m p Operation, and Treatment costs
Operational Cost
-
Jtin
(10)
i=1
c2Q~(h~ - zg8) -}- E c3Qi + Z(c5 - C6QT) dt i=ne +1
(11) where" n
=
total number of injection/extraction wells;
~e
=
number of extraction wells;
Qi
injection/extraction rate for well i (m3/s);
-
-
hmin Zgs
di tin, t f in
Z
Z
design extraction rate for well i (m3/s);
-
total extraction rate (m 3/s);
-
hydraulic head in well i (m); minimum allowable hydraulic head (m);
-
-
-
-
elevation of ground surface (m); depth of well i below the ground surface (m); initial and final time for remediation (s); air-stripping tower height (m);
]Q~[b5
[Ce--CElb6c?}
-
max{b3H b4
=
Henry's law coefficient = 0.2;
,
_> c E ) ;
influent concentration (kg/m3); CE
=
target effluent concentration = 5 10 .5 k g / m 3.
The coefficients bi and ci associated with the OF may be found in [8].
1167 4.2. C o n s t r a i n t s The state and decision variables are constrained by both hydraulic and solute mass constraints [8]. Hydraulic constraints include limitations on the extraction and injection rates at each well"
o ~e xa x <- -Q i < o- -m a xi n
i -- 1, 2, .., n
(12)
n
OT~ X
< - - ~-~Oi
(13)
i=1
In (12) ~ 0 (~max
max <
and -~i~ 0 "~x are the maximum extraction and injection rates at any well
0
--
" (~max ~
>
_
0
" [(~max
I __ Io m a x I - 6.410 "* in
-3
I113
/s
In (13) I O m a x is the maximum total extraction rate which has been set to -3 210 .2 m3/s Constraints on pumping rates (Equations (12)-(13)) are established by setting lower and upper bounds to the values of the decision variables between which the GA searches for the optimum. Hydraulic constraints include also limitations on the aquifer drawdown" '~T
~
hmin <__hi < hmax i -
"
"
(14)
1, 2, .., n
where h,~i~- 10 m and hmax=Zgs. Solute mass constraints are specified by requiring (a) the contaminant concentration to be zero on the fixed head boundaries shown in Figure 1 at any time during remediation; (b) a minimum aquifer cleanup performance as quantified by the global fraction of mass that may remain in the system at the end of the remediation horizon" C(x-
1000, y, z, t ~ < t <_ t f ~ ) -
C(x, y -
1000, z, t ~ < t <_ t f ~ ) -
0
(15) (16)
with the operation time t fin-tin equal to 10 years, and the PAT goal M set to 0.05. Both drawdown (Equation (14)) and solute mass constraints (Equations (15)-(16)) are enforced with high-order-of-magnitude penalty coefficients that simply multiply the total cost function if the constraints are not met. 5. R E S U L T S
AND
DISCUSSION
The management problem is solved by linking a GA code [10,11] to the F E / R W flow and transport model described in Section 2, which represents the "function call" calculating the cost associated with each particular design. In the equivalent-K (homogeneous) test case the optimum was found with about 470 OF evaluations. Figure 7 shows the optimal solution determined by the GA. Only 4 of the 8 potential wells are included in the optimal design. The maximum flow rate at well 8 is needed in order to prevent contaminant from reaching the upper boundary of the domain as prescribed by the constraint of Equation (15). Lower flow rates are needed from the other wells in the higher concentration areas of the plume so as to satisfy the minimum cleanup performance as in Equation (16).
1168
Figure 7. GA results for the 8-well design: Q2=Q6=4.63010-4m3/s (40 m3/day); Q4=2.31510-4m3/s (20 m3/day); Qs-3.47210-3m3/s (300 ran/day); wells 1, 3, 5, 7 are inactive.
Figure 8. Distribution of contaminant particles in the layered aquifer at the end of the 10-years long remedial time, where optimal policy obtained from equivalentK (homogeneous) case is applied to the layered case.
Table 2 reports the costs and the solute mass remaining associated with this solution. Table 2 also reports the results obtained when the optimal remediation policy determined with the equivalent-K (homogeneous) aquifer is applied to the heterogeneous aquifer. As anticipated, the performance of the policy calculated for the equivalent homogeneous aquifer is very poor. About 32% of the solute mass remains in the aquifer at the end of the remedial horizon. Figure 8 describes the distribution of contaminant particles in the layered aquifer 10 years after operations began. Solute remains primarily in the low permeable layer between the elevations of 15 and 18 In and the flow rates obtained from the equivalent-K case are not adequate.
Table 2 Comparison of remediation costs and performances. Cost ($) Aquifer Mr (%) Capital Pump Treat Homogeneous 0.5 182,380 6,790 63,770 Layered 32.6 185,630 6,960 65,740
Total 252,940 258,330(*)
(*) The cost does not reflect the p e n a l t y to be applied because constraint (16) is not met.
1169
Figure 9. Total mass of solute (kg) in the aquifer vs time. The period of time spans over 15 years from the start of contaminant release up until the end of the remedial operations.
In Figure 9, the behaviors of the total solute mass in the aquifer vs. time are displayed. The solid profile represents the plume development as simulated for the heterogeneous aquifer. The dashed profile indicates how contaminant is extracted from the aquifer following the optimal solution (Figure 7) determined assuming an equivalent homogeneous porous medium. The substantially double-slope shape of the diagram suggests that a solution with variable pumping rates would probably allow reducing the total cost. Finally, the dotted profile reveals how the contaminant would be cleaned up if the optimal policy of Figure 7 were applied to the layered heterogeneous aquifer. These results, though obtained using two simplistic conceptual models for the aquifer, the homogenous and the layered one, prove how important site characterization and collection of data may be in order to determine the optimal solution for the PAT process. We note, however, that the quantitative results presented here may be somewhat biased by neglecting the dispersion and soil moisture gradient terms in Equation (6).
1170 6. C O N C L U S I O N S
In this work we analysed how simplifying hypotheses may affect the final optimal remediation policy. A heterogeneous unconfined layered aquifer was compared to a homogeneous one with an equivalent (lumped) hydraulic conductivity. An optimal design was obtained for the homogeneous aquifer and this design was applied to the layered aquifer. Analysis of the two solutions, compared in terms of cleanup performance associated with the realization of the process, revealed that the performance of the policy calculated for the equivalent homogeneous aquifer results was in fact very poor. At the end of the planned remedial horizon, more than 30% of the solute is still in the aquifer, located in low permeable layers below the water table. The cleanup of these ~'stagnant" regions of the aquifer thus requires pumping rates and operation times much higher than those estimated in the homogeneous scenario. These results prove the importance of site characterization in order to determine the optimal solution for the PAT process. Future development of this work will consider the inclusion of site characterization in the objective function to evaluate when a more thorough knowledge of the aquifer is justified in term of both cost and remediation efficiency; the presence of both injection and extraction wells; the comparison of results for layered aquifers versus aquifers with higher degrees of heterogeneity; and the inclusion of dispersion and soil moisture gradient terms. REFERENCES
1. J . R . Philip, Adv. Hydrosci., 5 (1969) 215. 2. P.S. Huyakorn, S. D. Thomas, and B. M. Thompson, Water Resour. Res. 20 (1984) 1099. 3. P.S. Huyakorn and G. F. Pinder, Computational Methods in Subsurface Flow, 1983, Academic, San Diego. 4. C. Paniconi and M. Putti, Water Resour. Res. 30 (1994) 3357. 5. J. Bear, Hydraulics of Groundwater, 1979, McGraw-Hill, New York. 6. G. Uffink, NATO Advanced Workshop on Advances in Analytical and Numerical Groundwater Flow and Quality Modeling, Portugal, 1987. 7. A . F . B . Tompson, E. G. Vomvoris and L. W. Gelhar, Numerical Simulation of Solute Transport in Randomly Heterogeneous Porous Media: Motivation, Model development, and Application. Rep. UCID-21281, Lawrence Livermore Natl. Lab., Livermore, California, (1987). 8. A.S. Mayer, C. T. Kelley and C. T. Miller, Adv. Water Resour. 25 (2002) 1233. 9. E . R . LaBolle G. E. Fogg and A. F. B. Tompson, Water Resour. Res. 32 (1996) 583. 10. D. E. Golberg, Genetic algorithms in search, optimization, and machine learning, 1989, Addison-Wesley, Reading, MA. 11. M. Erickson, A. S. Mayer and J. Horn, Adv. Water Resour. 25 (2002) 51.