- Email: [email protected]
Vertical transport processes in unconfined aquifers
Journal of Contaminant Hydrology, 4 (1989) 93-107 Elsevier Science Publishers B.V., Amsterdam - - Printed in The Netherlands
93
VERTICAL TRANSPORT PROCESSES IN UNCONFINED AQUIFERS
DAVID W. OSTENDORF, DAVID A. RECKHOW and DAVID J. POPIELARCZYK Environmental Engineering Program, Civil Engineering Department, University of Massachusetts, Amherst, MA 01003, U.S.A. (Received August 18, 1987; revised and accepted April 12, 1988)
ABSTRACT Ostendorf, D.W., Reckhow, D.H. and Popielarczyk, D.J., 1989. Vertical transport processes in unconfined aquifers. J. Contam. Hydrol., 4: 93-107. We derive simple two-dimensional mathematical models describing the unsteady transport of conservative contaminants through an unconfined aquifer with a gently sloping aquiclude subject to advection, recharge, and vertical dispersion. The inclusion of vertical transport terms permits the proper nonreactive analysis of closed and open chemical systems, with the latter allowing dispersion of volatile constituents across the water table. These systems exhibit conservative and pseudoreactive behavior respectively when the pollution is analyzed on a depth-integrated basis, as is common in present one-dimensional models of groundwater contamination. Vertical and longitudinal chloride and total inorganic carbon observations at the well-documented Babylon, Long Island sanitary landfill plume are used to calibrate and test the analyses with a modest level of accuracy, using the vertical dispersivity as a calibration factor in this testing process. The parameter is important in the determination of reaeration rates across the water table and nutrient mixing from below in the related problem of biological transformations near the free surface. INTRODUCTION W e m o d e l t h e t w o - d i m e n s i o n a l u n s t e a d y t r a n s p o r t o f a n o n r e a c t i v e constituent through an unconfined aquifer with a gently sloping planar bottom. The analysis proceeds under steady hydraulics and includes an approximate account of the coupled vertical transport processes of recharge and dispersion, subject to closed (nondispersive) and open (dispersive) boundary constraints at the water table. The latter case governs the loss of volatile constituents through the free surface. The resulting decrease of concentration appears as a reactive loss term, or sink, in conventional depth-integrated models of s u b s u r f a c e p o l l u t i o n ( W i l s o n a n d M i l l e r , 1978; B e a r , 1979; P r a k a s h , 1982). Hence the new analysis distinguishes the non-reactive degassing of such species as carbon dioxide or nitrogen from truly reactive phenomena in which the losses are intrinsic to the chemistry and microbiology of the subsurface e n v i r o n m e n t . T h e r e c h a r g e a n d d i s p e r s i v i t y a r e f u n c t i o n s o f t h e flow field, a n d accordingly may be used to characterize the related and important processes of aquifer reaeration and nutrient mixing with depth below the water table. These
0169-7722/89/$03.50
© 1989 Elsevier Science Publishers B.V.
94 latter phenomena are of interest to investigators of subsurface microbiology and in situ hazardous waste treatment. Our model builds upon the existing analytical modeling efforts of Gelhar and Wilson (1974) and Ostendorf et al. (1984), who derive the recharge hydraulics and source term behavior for the present study, as discussed below. The advectire-dispersion solutions of Ogata and Banks (1961) and Parker and Van Genuchten (1984) describe the open- and closed-system profiles respectively, although the transport is horizontal in the original contexts. Additional, less directly related, depth-integrated models may be cited as well: Lenau (1972) postulates a steady-state, conservative injection of pollution from a recharge well in a confined aquifer, while Wilson and Miller (1978) consider the transient response of a simply reactive contaminant from a similar source configuration. Chen (1987) derives a model including radial hydraulics and transport of a conservative species in a confined aquifer. Bear (1979) summarizes unsteady contaminant migration due to a series of one-dimensional reactive sources, and Prakash (1982) models steady-state reactive pollution due to instantaneous point, line, and volume sources. The method of characteristics adopted here follows one-dimensional analytical applications in the unsaturated and saturated zones by Wilson and Gelhar (1981) and Charbeneau (1981), respectively. Valocchi (1985, 1986) examines the validity of the local equilibrium assumption for simply adsorbing solutes in one-dimensional and radial flow fields on an analytical basis as well. The present simple analytical modeling approach complements existing more complex numerical contaminant transport codes. In this regard, we may distinguish two (Konikow and Bredehoeft, 1978; Borden and Bedient, 1986) and three (Huyakorn et al., 1986) -dimensional numerical solute transoirt models of single and coupled (Cederberg et al., 1985) constituents. These site-specific programs, with their attendant documentation requirements, are appropriate for the detailed analysis of fully measured plumes in complex geological settings. As such, they contrast with the present analytical effort, which is intended to provide a readily used, generic account of two-dimensional contamination of an idealized or sparsely measured flow field by a single species. HYDRAULICS We begin by considering the two-dimensional conservation of water mass in the unconfined aquifer, wherein a small (near)-vertical, average linear velocity component w is imposed by recharge upon a vertically uniform, (near)horizontal, average linear velocity u (Gelhar and Wilson, 1974): ~u ~x
--
+
0w ~z
-
0
(1)
with distance z below the water table and downgradient distance x from the source of pollution, as sketched in Fig. 1. Following Ostendorf et al. (1984), the unconfined aquifer flow field is perturbed by small changes due to a gently
95 landfill
x x x
x×x
It ~,w hs L _ _ ] -
~
source plane
----,
I
]
x,u h
Fig. 1. Definition sketch. sloping aquiclude, head loss, and recharge:
[ X(n
u = us 1 + ~ tanfl = tanfl'
s
-
tanfl
)1
nvus
kg
(2a) (2b)
with aquifer depth h, kinematic viscosity v, permeability k, gravitational acceleration g, recharge e, and porosity n. The s subscript denotes pollutant source conditions at the downgradient boundary of the landfill, while the slope tan fl relative to the water table includes the horizontally based slope tan fl'. The corresponding vertical average linear velocity follows from eqns. (1) and (2a). To leading order, we have: w -
n
- us tan
hs
(3)
Equation (3) suggests that w equals e / n at the water table and yields a velocity parallel to the aquielude where z equals ha, the first-order estimate of the actual aquifer thickness h. We are interested in frames of reference moving at speed u from the downgradient end of the landfill, taken as the (far field) source of pollution, as indicated by Fig. 1. The velocity and path xf of the frames are given by Ostendorf et al. (1984) as: dxf dt
--
(4a)
U
xf[
--tan )l
xf(~
z = -- 1 -U~ ~ss
rtUs
(4b)
96 =
t - ts
(4c)
with time r in the m o v i n g r e f e r e n c e frame which leaves the s o u r c e at time t~. A (near field) l i n e a r r e s e r v o i r analysis u n d e r the landfill has been derived by O s t e n d o r f et al. (1984) in o r d e r to specify s o u r c e plane c o n c e n t r a t i o n as a f u n c t i o n of ts and user population. Following the c o n c e p t of G e l h a r and Wilson (1974), this simple a p p r o a c h t r e a t s the complex n e a r field flow region as a fully mixed c h e m o s t a t r o u t i n g pollution input from the landfill to the o u t p u t plane at the d o w n g r a d i e n t end of the facility, w h e r e it serves as a v e r t i c a l l y uniform source of c o n t a m i n a t i o n . Thus, the frames c a r r y k n o w n pollution levels into the far field with presumed uniform initial v e r t i c a l profiles. OPEN- AND CLOSED-SYSTEM CONCENTRATION PROFILES The c o n s e r v a t i o n of n o n r e a c t i v e c o n t a m i n a n t mass c o n c e n t r a t i o n c can be c o n v e n i e n t l y studied in the m o v i n g r e f e r e n c e frame as a b a l a n c e of storage, v e r t i c a l advection, and v e r t i c a l dispersion: e~c ~ ~c ~2c ~-~ + -n--~z - D--0z2 = 0
(5a)
D
(5b)
= ~us
L o n g i t u d i n a l dispersion is neglected in this t r a n s p o r t e q u a t i o n on the premise of an appreciable l o n g i t u d i n a l velocity and a c o n t i n u o u s pollution source ( O s t e n d o r f et al., 1984); ~ a c c o r d i n g l y r e p r e s e n t s the v e r t i c a l dispersivity (Bear, 1979). The use of sin in place of w in the advective t r a n s p o r t term reflects o u r i n t e r e s t in c o n d i t i o n s n e a r the w a t e r table [see eq. (3)], w h e r e the bulk of t r a n s p o r t o c c u r s in the form of dilution and degassing t h r o u g h the r e c h a r g e lens. The g o v e r n i n g c o n s e r v a t i o n e q u a t i o n (5a) is solved subject to a presumedly uniform s o u r c e p l a n e c o n c e n t r a t i o n cs and a deep aquifer, i.e.: c =
cs
(T = 0)
(6a)
c =
cs
(z =
(6b)
~)
The l a t t e r c o n s t r a i n t suggests t h a t most of the t r a n s p o r t occurs n e a r the free surface and is u n a f f e c t e d by the bottom; an assumption in keeping with the use of E/n as the v e r t i c a l v e l o c i t y estimate. The r e m a i n i n g w a t e r table b o u n d a r y c o n d i t i o n differs for open and closed chemical systems: c =
ca
0c ~ - Dc~--~ + -nC = -nCa
(z = 0, open)
(7a)
(z = 0, closed)
(7b)
with c o n s t a n t a m b i e n t c o n c e n t r a t i o n c~. T h e open-system c o n d i t i o n (7a) rests on the premise of a rapidly diffusing u n s a t u r a t e d zone, as discussed for i n o r g a n i c c a r b o n in the Appendix.
97 Equations (5a), (6) and (7a) constitute the classical advective-diffusion problem with solution (Ogata and Banks, 1961): c~
ca
2
erfc L 2 ~ - - ~ J
+ exp ~-~ erfc L 2(Dr),. 2 j j
(8a)
(open) erfc y -
it,/2
exp ( - y ' 2 ) d y '
(8b)
Y
with complementary error function erfc y (Abramowitz and Stegun, 1972). Asymptotic behavior of this function may be used to advantage, whence: (z - ~ / n ) 2 7
c - Ca -- 1 -- i { rz__ ~,/n 1 (_~y,2exp ca - ca ~ erfcL2(Dr)l/2_j + 2
z + ~r/n
-4D-~T] ]1 z + ~ j
> 5(Dr) '/2
(9a) (9b)
The closed-system condition (7b) is predicated on a no-flux constraint at the water table, so that the advective and dispersive mechanisms are in balance at that location. A more complicated solution follows, as derived by Parker and Van Genuchten (1984):
c- Ca-
] + exp(~'Z )
cs -- ca
[Z +---~.1[~
~
( ~ ~1,,~ n \~)
[
erfc L 2(Dr)1/2 J
(z - ~/n) ~]
exp
]~-
+
1
~(z + 8"In)1 2On ] [~_ : ~/n l
-j - ~ erfc L 2(Dr) '/2 J
(10)
This too simplifies for a large argument: (z - ~r/n)'-']
C--Ca --
1+
c~ - ca
z + ~r/n
JlerfcrZ-
exp
z + ~:~/n > 5(Dr) 1/2
2
/n 1
(lea)
L2(-D~-~J (llb)
Figure 2 displays typical profiles for open and closed systems under dominant dispersion and recharge conditions for a pure ambient aquifer. The strong dispersion case features deeper mixing of open constituents and an appreciable free-surface closed-system concentration. Strong recharge on the other hand, establishes a pure overlying lens of water regardless of the open or closed nature of the pollutant. The area under the contaminant profiles is proportional to the depth averaged concentration: Fig. 2a indicates that open and closed systems exhibit values below their source levels. The open-system decrease is biggest since it includes dilution and degassing effects. The dilution reduction is predominant in Fig. 2b and is common to both systems.
98 0
~"--" ope n
-,~sed
0
"~, ~
\ open and closed
Z
hs
Z
0.5
hs 0.5
I
(°)
I.O o
c
o.5 / cs
,.o (b)
I
0
0,5
I.O
1.0
c~ cs
Fig. 2. Typical open and closed c o n t a m i n a n t system profiles for vs = 10-Sms -1, n = 0.3, and = 10Ss. S t r o n g dispersion (e = 1 0 - g m s -1, a = 10-~m) is sketched in Fig. 2a, while s t r o n g r e c h a r g e (~ = 10-Sins -~, ~ = 10-3m) is s h o w n in Fig. 2b. CALIBRATION AND T E S T I N G
The total inorganic carbon and chloride observations of Kimmel and Braids (1974, 1980), supplemented by raw data supplied by the USGS may be used to calibrate and test the open and closed models, respectively. Vertical and longitudinal profiles of contamination from a sanitary landfill on Babylon, Long Island provide ideal and independent data bases, since the hydraulics are well understood and the unconfined aquifer geometry is appropriately simple. Indeed, the existing depth-integrated model of Ostendorf et al. (1984) is calibrated with this plume, based on a source average linear velocity (us) of 3.37 × 10 6ms-l, recharge (e) of 3.25 × 10 9ms 1, aquifer depth (ha) of 22.5m at the source plane, permeability of 6.34 × 10-11m ~, viscosity of 1.1 × 10-6m2s -1, and a porosity of 0.27. Source plane concentrations are computed from monitoring well data located at the downgradient boundary of the landfill, following the near field procedure proposed by Ostendorf et al. (1984). The source plane wells are used to calibrate a per capita contaminant generation rate Scl equals 1.40 × 10 _8 kg c a p - i s -1 for chloride in the earlier analysis, and this value is adopted here to estimate the ca values in closed system testing. A total inorganic carbon value of Sc equals 7.90 × 10-gkgcap-ls -1 generates the open-system source term concentrations in the present study, based on 1973 observations in the source plane monitoring wells. The new Sc value is needed since Ostendorf et al. (1984) studied bicarbonate and not total inorganic carbon in their prior work. We compare data and theory using statistics of the error 5 defined by: 5 =
c (measured) - c (predicted) c (measured)
(12a)
= _1.Z5 ]
(12b)
= ~ . Z&2 - ~ t I/2
(12c)
99
T h e m e a n e r r o r ~ and s t a n d a r d d e v i a t i o n a are c o m p u t e d in a c c o r d a n c e with B e n j a m i n and Cornell (1970). The sign of 3 indicates model over or underprediction and is a c c o r d i n g l y useful in identifying systematic model errors. The e r r o r s t a n d a r d d e v i a t i o n is based on the absolute v a l u e of individual 3's and conseq u e n t l y m e a s u r e s the m a g n i t u d e of the error. In this regard, a b o u t 2/3 of our predictions lie w i t h i n a of t h e i r m e a s u r e d values for a zero m e a n error. T h e (far field) l o n g i t u d i n a l total i n o r g a n i c c a r b o n profile d a t a observed in 1973 by Kimmel and Braids (1980) are used to c a l i b r a t e the vertical dispersivity in eqns. (5b) and (8a) with the results summarized in Table 1 and Fig. 3a. The TABLE 1 L o n g i t u d i n a l profile of total inorganic carbon: model calibration x (m)
~ (10 s s)
cs (kg m - 3)
z (m)
c (meas) (kg m 3)
c (pred) (kg m - 3)
(%)
360
1.05
0.136
630 900
1.80 2.54
0.117 0.091
920
2.59
0.089
1550 1890 2180
3.89 5.02 5.69
0.057 0.045 0.038
2230 2260 2810
5.80 5.87 7.05
0.037 0.036 0.020
1.75 9.05 25.11 23.16 3.08 8.97 18.45 11.74 21.50 12.17 12.69 8.75 19.33 17.66 21.38 11.89 18.26
0.017 0.087 0.141 0.091 0.007 0.136 0.137 0.035 0.060 0.055 0.067 0.018 0.037 0.010 0.040 0.010 0.010
0.022 0.100 0.136 0.114 0.015 0.044 0.077 0.055 0.081 0.029 0.020 0.011 0.024 0.021 0.024 0.007 0.010
- 32 - 15 4 24 - 115 68 44 56 - 36 47 69 42 36 99 40 32 - 6
Notes: 1973 data (Kimmel and Braids, 1980), supplemented by USGS r a w data. Source conditions from O s t e n d o r f et al. (1984). C o n c e n t r a t i o n s expressed as carbon. Calibrated dispersivity ~ 0.073 m, zeros m e a n error. Calibration e r r o r s t a n d a r d deviation a = 54%. 0.15
el
I
0.10 "u'x"~. ~
c (kg m-3)
I I0
"a = 0.073 m
Z (rn)
o
0.05__
0
I
©15-25m o5-15m eO-5m
(0) 0
II
20
(b) I xaXt
50 500
1000 x
1500 (m)
2000
2500
3000
0
l
0.025 0.050 0.0"75 ¢(kg m -3)
0.0[0
Fig. 3. Total c a r b o n calibration (a) and test (b) data (circles) and predictions (lines). 1973 observations of Kimmel and Braids (1974, 1980), supplemented by raw data supplied by USGS.
100
more distant wells have larger travel times (~), corresponding to frames leaving the landfill at earlier times with smaller user populations hence lower source concentrations. The vertical location refers to distance below the average water table elevation, and the open system is chosen for calibration purposes since it possesses a steeper gradient than its closed count erpart and is accordingly more sensitive to ~. The calibration is a Fibonacci search (Beveridge and Schechter, 1970) for the dispersivity value zeroing the mean error. The operation yields a standard deviation of 54% associated with the optimal dispersivity: = 0.073 m
(13)
This dispersivity value is close to the first-order decay based value of 0.10 m put forth by Ostendorf et al. (1984) on an ad hoc basis in their depth-integrated analysis of bicarbonate data at the site. The newer version reflects a more rigorous estimate of the transport mechanisms: the "decay" is in reality a mass transport process constrained by vertical mixing and recharge. Table 1 suggests no systematic dependence of the error upon horizontal or vertical distance, so that the relatively high standard deviation may in part be attributed to sampling errors. In any event, additional independent testing of the calibrated model is certainly warranted. The vertical inorganic carbon, longitudinal chloride, and vertical chloride observations of Kimmel and Braids (1974, 1980) provide these necessary tests, as summarized in Tables 2-4 and Figs. 3b and 4. The accuracy of the vertical inorganic carbon test is encouraging, as evidenced by the 28% mean and 21% standard deviation error statistics. The longitudinal chloride data provide an additional measure of support, particularly if one neglects two data points immediately beneath the water table. Most of the mean error ( - 4 6 % ) and standard deviation (126%) for this test is created by the two near-surface data points 360m and 900m downgradient: we compute quite accurate values of - 2°./o and 25% for the statistics without these samples. The vertical chloride TABLE 2 Vertical profile of total inorganic carbon: model test z (m)
c (measured) ( k g m ~)
c (predicted) ( k g m -3)
3 (%)
5.58 12.01 14.27 19.30 24.17
0.0152 0.0382 0.0405 0.0921 0.0420
0.0113 0.0249 0.0294 0.0378 0.0442
25 34 27 59 - 5
Notes: USGS r a w data, except 24.17 m value (Kimmel and Braids, 1974). Source condition from O s t e n d o r f et al. (1984). C o n c e n t r a t i o n s expressed as carbon. Profile location x = 1570m, r = 4.26 × 108s. Test statistics ~ = 28%, a = 21%.
101 TABLE 3 L o n g i t u d i n a l profile of chloride: model t e s t x (m)
~ (10Ss)
c~ ( k g m 3)
z (m)
c (meas) ( k g m -3)
c (pred) ( k g m 3)
(%)
360
1.05
0.258
630 900
1.80 2.54
0.224 0.177
920
2.59
0.174
1550 1890 2180
3.89 5.02 5.69
0.117 0.096 0.087
2230 2260 2810
5.80 5.87 7.05
0.080 0.079 0.051
1.75 9.05 25.11 23.16 3.08 8.97 18.45 11.74 21.50 12.17 12.69 8.75 19.33 17.66 21.38 11.89 18.26
0.054 0.190 0.270 0.180 0.024 0.200 0.360 0.130 0.180 0.109 0.150 0.060 0.057 0.057 0.067 0.042 0.045
0.213 0,247 0.258 0.223 0.134 0.156 0.172 0.159 0.171 0.102 0.081 0.067 0.078 0.071 0.072 0.042 0.045
- 294* - 30 5 - 24 - 459* 22 52 - 22 5 6 46 - 12 - 38 - 25 - 8 0 0
Notes: 1973 d a t a ( K i m m e l a n d Braids, 1980), s u p p l e m e n t e d by U S G S r a w data. S o u r c e c o n d i t i o n s from O s t e n d o r f et al. (1984). A m b i e n t c o n c e n t r a t i o n of 0.020 k g / m 3 included. T e s t s t a t i s t i c s w i t h n e a r - s u r f a c e p o i n t s ~ = - 46%, ~ = 126%. * T e s t s t a t i s t i c s w i t h o u t n e a r - s u r f a c e p o i n t s ~ = - 2%, ~ = 25%.
TABLE 4 V e r t i c a l c h l o r i d e profile: model t e s t z (m)
c (measured) ( k g m -3)
c (predicted) ( k g m -3)
(%)
5.58 12.01 14.27 19.30 24.17
0.088 0.140 0.170 0.200 0.170
0.082 0.094 0.097 0.102 0.106
6 33 43 49 39
Notes: U S G S r a w data, e x c e p t 24.17m d a t a ( K i m m e l a n d Braids, 1974). S o u r c e c o n d i t i o n s from O s t e n d o r f et al. (1984). Profile l o c a t i o n x = 1570m, ~ = 4.26 × 108s. T e s t s t a t i s t i c s ~ ~ 34%, ~ = 15%. test summarized the
34%
ambient
mean chloride
the plume.
in Table and
15%
4 and Fig. 4b does not fare as well however, standard
concentration
The closed-system
deviation.
of 0.020 kgm errors
do suggest
Both
chloride
-3, as indicated a tendency
tests
based
on
reflect
an
by wells
to overpredict
outside near
102 0.40
I
01
I
I
,,5!25m o5-15m •O-5m
0.30 c (kg m-3)
o % .
0.20
0
I
ct= 0.073 m
%.
- -
I0
Z (m)
2O
0. I0
(o)
I 500
0
"i I000
1
I
1500
2000
I
(b)
30
2500
0
3000
I
0.05
I
O.lO 0.15 c (kg m-3)
0.20
X (m]
Fig. 4. Longitudinal (a) and vertical (b) chloride test data (circles) and predictions (lines). 1973 observations of Kimmel and Braids (1974, 1980), supplemented by raw data supplied by USGS. the free surface and u n d e r p r e d i c t at d e e p e r locations in the profile, perhaps in response to a n o n u n i f o r m s o u r c e b e h a v i o r not included in the present analysis. T a k e n as a whole then, the c a l i b r a t i o n and test results provide a modest level of s u p p o r t for this analysis, a l t h o u g h a d e p t h - v a r y i n g s o u r c e term c h a r a c t e r i z a tion could well improve the model a c c u r a c y for the closed-system constituents. RESULTS T h e c a l i b r a t e d and tested open-system model m a y be c o m p a r e d with its one-dimensional c o u n t e r p a r t to e l u c i d a t e the physical basis of the ad hoc decay c o n s t a n t A used in c o n v e n t i o n a l models. T h e d e p t h - a v e r a g e d c o n c e n t r a t i o n follows upon i n t e g r a t i o n of eqn. (8a), as derived in the Appendix: ca cs - c.
1
ierfc y
~
=
h,
ierfc
erfc(y')dy"
':2]
-
erf
k
]
(14a)
(14b)
Y
erfy
=
1 - ierfcy
(14c)
with e r r o r f u n c t i o n i n t e g r a l ierfc y defined by A b r a m o w i t z and S t e g u n (1972). E q u a t i o n (14) m a y be used with d e p t h - i n t e g r a t e d d a t a to c a l i b r a t e a vertical dispersivity v a l u e in the absence of m o r e definitive v e r t i c a l open-system conc e n t r a t i o n profiles. T h e " t r u e " one-dimensional model r e p r e s e n t e d by eqn. (14) c o n t r a s t s with the c o n v e n t i o n a l b a l a n c e of a d v e c t i o n and first-order d e c a y (Bear, 1979): - ca =
(cs - c , ) e x p ( - / l ~ )
(15)
in t h a t r e c h a r g e and v e r t i c a l dispersivity p a r a m e t e r s a p p e a r explicitly in the f o r m e r e x p r e s s i o n as a c o n s e q u e n c e of the p r i o r d e p t h - v a r y i n g analysis. The ad
103
hoc decay constant can be related directly to the responsible vertical transport properties for small ~ values, whence: --
Ca
Cs - -
Ca
(1 - ~tr)
- C, cs - Ca
1 -
2 I__~l I/2 ~
(ad hoc decay)
(16a)
(true mixing)
(16b)
Thus dispersion induces pseudoreactive decay of a degassing constituent, with an ad hoc constant approximately given by: -
)+ ~
(17a)
-
h+ \ ~ X c /
xc = u+~
(17b)
with the plume length scale Xc in the longitudinal direction specified by advection in the unconfined aquifer. We may also estimate the flux J of open system contaminant through the water table once the concentration profile is specified, since J is the combined result of vertical recharge advection and dispersion (Bear, 1979): J
=
nw(c-
c.) -
Oc
nD-~z
(z = 0)
(18)
In view of eqns. (7a) and (8), we find:
g = -(c+-
{i/D'~l/2 [- (~.c/n)2_l
ca) n[~)__
exp|_
4D+ J -
~
-/"
~.1/2
x)
~erfc[~)~
(19)
The contrasting roles of the vertical transport mechanisms may be noted for limiting cases of strong recharge and strong dispersion J
~
-n(c~
J
~ 0
(O~']2
- c . ) \-~-~/
=
~)
(20a)
(s =
~)
(20b)
(D
Dispersion mixes the carbon towards the water table at a rate that decreases with travel time due to a flattening vertical gradient. The vertical recharge flow flushes pollution back into the plume by downward advection, and so reduces J. CONCLUSIONS
We derive simple two-dimensional mathematical models describing the unsteady transport of conservative contaminants through an unconfined aquifer with a gently sloping aquiclude subject to advection, recharge, and vertical dispersion. Open (dispersive) and closed (nondispersive) chemical
104 systems are considered by imposition of appropriate boundary conditions at the water table, and the pseudoreactive depth-integrated behavior exhibited by the volatile open-system constituents is explained in terms of fundamental vertical mixing parameters. Total inorganic carbon and chloride data at the Babylon, Long Island landfill are used to calibrate and test the open- and closed-system models with a modest level of accuracy. The observed longitudinal distribution of total inorganic carbon is used for calibration, an exercise yielding a vertical dispersivity estimate of 0.073m with a standard deviation of 54%. Vertical profiles of inorganic carbon and chloride as well as a longitudinal chloride data set are used to test the theory, an exercise resulting in means and standard deviations ranging from - 2 to 34% in magnitude. The statistics do not include two chloride data points immediately below the water table, which are substantially overpredicted by the model. Future research may proceed along various related fronts. The present analysis of groundwater contamination rests on a sequential understanding of the following phenomena, both in terms of data and theory: - - groundwater hydraulics - - pollution source history depth-integrated contamination, conservative contaminants depth-varying contamination, conservative contaminants The first extension of the modeling effort should be the inclusion of a depthvarying source condition in response to the systematic errors noted in chloride testing. Such a modification should preserve the simplicity inherent in an analytical model approach while maintaining the accuracy of the open-system profile tests in the present study. Additional extensions of the modeling effort might include the depth-varying transport of truly reactive pollutants, first for a single species and then for coupled contaminants. One is certainly drawn to dissolved oxygen as a reactive constituent of interest, and the speciation of inorganic carbon at Babylon seems a logical starting point for the coupled work. The present model should be modified to account for the influence of a nondispersive bottom boundary to expand its applicability, with the possible inclusion of an appropriately simple allowance for the vertical variation of the vertical velocity w. Additional testing of the present theory is certainly called for, at other sites and for other contaminants. In the latter regard, total nitrogen in a sewage effluent plume seems a likely candidate for open system testing, due to degassing associated with the denitrification process. -
-
-
-
ACKNOWLEDGMENTS Partial funding for this research was provided by the Massachusetts Division of Water Pollution Control under Grant Number 86-30765; the Authors acknowledge and appreciate their support. We also thank the Water Resources Division of the United States Geological Survey for providing the raw data base for the Babylon landfill.
105 REFERENCES Abramowitz, M. and Stegun, I.A., 1972. Handbook of Mathematical Functions. National Bureau of Standards, Washington, DC, 1046 pp. Bear, J., 1979. Hydraulics of Groundwater. Wiley, New York, NY, 567 pp. Benjamin, J.R. and Cornell, C.A., 1970. Probability, Statistics, and Decision for Civil Engineers. McGraw-Hill, New York, NY, 684 pp. Beveridge, G.S.G. and Schechter, R.S., 1970. Optimization: Theory and Practice. McGraw-Hill, New York, NY, 773 pp. Borden, R.C. and Bedient, P.B., 1986. Transport of dissolved hydrocarbons influenced by oxygen limited biodegradation 1. Theoretical development. Water Resour. Res., 22: 1973-1982. Cederberg, G.A., Street, R.L. and Leckie, J.O., 1985. A groundwater mass transport and equilibrium chemistry model for multicomponent systems. Water Resour. Res., 21:1095 1104. Charbeneau, R.L., 1981. Groundwater contaminant transport with adsorption and ion exchange chemistry: method of characteristics for the case without dispersion. Water Resour. Res., 17: 705 713. Chen, C.S., 1987. Analytical solutions for radial dispersion with Cauchy boundary at injection well. Water Resour. Res., 23:1217 1224. Gelhar, L.W. and Wilson, J.L., 1974. Groundwater quality modeling. Ground Water, 12: 399~408. Hillel, D., 1982. Introduction to Soil Physics. Academic Press, New York, NY, 364pp. Huyakorn, P.S., Jones, B.G. and Andersen, P.F., 1986. Finite element algorithms for simulating three-dimensional groundwater flow and solute transport in multilayer systems. Water Resour. Res., 22:361 374. Kimmel, G.E. and Braids, O.C., 1974. Leachate plumes in a highly permeable aquifer. Ground Water, 12: 388-393. Kimmel, G.E. and Braids, O.C., 1980. Leachate plumes in groundwater from Babylon and Islip landfills, Long Island, New York. U.S. Geol. Surv., Prof. Pap. 1085, 38 pp. Konikow, L.F. and Bredehoeft, J.D., 1978. Computer model of two-dimensional solute transport and dispersion in groundwater. Techn. Water Resour. Invest. USGS, USGS, Reston, VA. Lenau, C,W., 1972. Dispersion from recharge well. J. Eng. Mech. Div., ASCE, 98: 331-334. Ogata, A. and Banks, R.B., 1961. A solution of the differential equation of longitudinal dispersion in porous media. U.S. Geol. Surv., Prof. Pap. 411-A, 7 pp. Ostendorf, D.W., Noss, R.R. and Lederer, D.O., 1984. Landfill leachate migration through shallow unconfined aquifers. Water Resour. Res., 20: 291-296. Parker, J.C. and Van Genuchten, M.T., 1984. Flux averaged and volume averaged concentrations in continuum approaches to solute transport. Water Resour. Res., 20: 866-872. Prakash, A., 1982. Groundwater contamination due to vanishing and finite size continuous sources. J. Hydraul. Div., ASCE, 108: 572-590. Stumm, W. and Morgan, J.J., 1981. Aquatic Chemistry. Wiley, New York, NY, 780 pp. Valocchi, A.J., 1985. Validity of the local equilibrium assumption for modeling sorbing solute transport through homogeneous soils. Water Resour. Res., 21: 808-820. Valocchi, A.J., 1986. Effect of radial flow on deviations from local equilibrium during sorbing solute transport through homogeneous soils. Water Resour. Res., 22: 1693-1701. Wilson, J.L. and Gelhar, L.W., 1981. Analysis of longitudinal dispersion in unsaturated flow 1. The analytical method. Water Resour. Res., 17:122 130. Wilson, J.L. and Miller, P.J., 1978. Two-dimensional plume in uniform groundwater flow. J. Hydraul. Div., ASCE, 104: 503-514. APPENDICES
Inorganic carbon boundary condition T h e i n o r g a n i c c a r b o n c o n c e n t r a t i o n c o n s i s t s o f b i c a r b o n a t e cl a n d d i s s o l v e d c a r b o n d i o x i d e c2 c o n t r i b u t i o n s ( S t u m m a n d M o r g a n , 1981) i n t h e p H r a n g e
106
(5-7) common to sanitary landfill plumes like t h a t at Babylon, Long Island: c = cl + c2
(21a)
cl = c2 ~
(21b)
K~ = 10-63mol L -1
(21c)
The equilibrium equation (21b) includes the hydrogen ion activity [H +] and equilibrium constant K1 expressed in moles/liter. The latter parameter is temperature and concentration dependent to a modest degree, but our intent here is an order of magnitude discussion of mass fluxes and a representative value is sufficient. The inorganic carbon and its species are expressed in terms of k g m -3 carbon, for consistency with the main analysis. The carbon dioxide species governs the total inorganic carbon degassing rate, since it is in gaseous equilibrium with the unsaturated zone air monolayer at the water table. Henry's law is appropriate (Stumm and Morgan, 1981) at this location, and may be expressed as:
c2 = KHRTp
(z = 0)
(22a)
R = 0.0821 atm L m o l - ' K 1
(22b)
KH = 0.045 mol L - ' a t m - '
(22c)
with Henry's constant KH, universal gas constant R, temperature T, and gaseous carbon dioxide density p, in k g m -3 carbon. For typical groundwater temperatures of 290 K, eqn. (22) suggests that the dissolved concentration and gaseous density are of equal magnitude on either side of the free surface: c2 ~
p
(z = 0)
(23)
The carbon dioxide (hence total inorganic carbon) flux J through the unsaturated zone is a gaseous dispersion process which may be estimated in accordance with (Hillel, 1982): J
~
-D~ P
(24)
with unsaturated zone thickness b and gaseous dispersivity Da of order 10- 6 m2s-, in magnitude. This crude estimate may be compared with its liquid side counterpart: J
~
c
-D-h~
(25)
in order to assess the magnitude of c at the water table. Collecting eqns. (21)-(25), we deduce: c(z = 0) c
D(
D~ 1 + ~
K1)
(26a)
107 C(Z = 0) ~
C
(pH
< 7.3)
(26b)
where the s a t u r a t e d and u n s a t u r a t e d zone thicknesses are the same order of m a g n i t u d e in size. P h y s i c a l l y speaking, the u n s a t u r a t e d zone efficiently carries off w h a t e v e r total carbon the s a t u r a t e d zone can deliver due to the smallness (1/10) of the dispersivity ratio; the simple b o u n d a r y condition (7a) is the result.
Open system The depth-integrated open-system c o n c e n t r a t i o n is defined as: hs
1 f cdz
-
(27)
hs
so that, recalling eqn. (8a): = c~ I =
i
cs - ca 2hs (I + II)
(28a)
L[zb/--n~--~/dz - ~ln]
(28b)
erfc
0 II --
i (sn)~ exp
L[z+~z/nl-2~D~-~ J
erfc . - d z
(28c)
0
The use of the infinite upper limit of i n t e g r a t i o n in eqn. (28) reflects the assumption of an infinite bottom b o u n d a r y condition in eqn. (6b). The first integral I is straightforward, while the second II requires an i n t e g r a t i o n by parts: II
nD~erfc
-
II -
+
nD [ erfc ( n i p2=~''cli2"~ ~, ]i ~ . -
(D~ 'i2
\~--~]
erfc
i[ exp
-_ aUn)~.~ 4D~
\2nD'121#_] (e'rli2"~]
j dz
(29a)
(29b)
This last relation simplifies with the aid of an identity: erfc(-y)
= 2 - erfcy
giving rise to eqn. (14a).
(30)
93
VERTICAL TRANSPORT PROCESSES IN UNCONFINED AQUIFERS
DAVID W. OSTENDORF, DAVID A. RECKHOW and DAVID J. POPIELARCZYK Environmental Engineering Program, Civil Engineering Department, University of Massachusetts, Amherst, MA 01003, U.S.A. (Received August 18, 1987; revised and accepted April 12, 1988)
ABSTRACT Ostendorf, D.W., Reckhow, D.H. and Popielarczyk, D.J., 1989. Vertical transport processes in unconfined aquifers. J. Contam. Hydrol., 4: 93-107. We derive simple two-dimensional mathematical models describing the unsteady transport of conservative contaminants through an unconfined aquifer with a gently sloping aquiclude subject to advection, recharge, and vertical dispersion. The inclusion of vertical transport terms permits the proper nonreactive analysis of closed and open chemical systems, with the latter allowing dispersion of volatile constituents across the water table. These systems exhibit conservative and pseudoreactive behavior respectively when the pollution is analyzed on a depth-integrated basis, as is common in present one-dimensional models of groundwater contamination. Vertical and longitudinal chloride and total inorganic carbon observations at the well-documented Babylon, Long Island sanitary landfill plume are used to calibrate and test the analyses with a modest level of accuracy, using the vertical dispersivity as a calibration factor in this testing process. The parameter is important in the determination of reaeration rates across the water table and nutrient mixing from below in the related problem of biological transformations near the free surface. INTRODUCTION W e m o d e l t h e t w o - d i m e n s i o n a l u n s t e a d y t r a n s p o r t o f a n o n r e a c t i v e constituent through an unconfined aquifer with a gently sloping planar bottom. The analysis proceeds under steady hydraulics and includes an approximate account of the coupled vertical transport processes of recharge and dispersion, subject to closed (nondispersive) and open (dispersive) boundary constraints at the water table. The latter case governs the loss of volatile constituents through the free surface. The resulting decrease of concentration appears as a reactive loss term, or sink, in conventional depth-integrated models of s u b s u r f a c e p o l l u t i o n ( W i l s o n a n d M i l l e r , 1978; B e a r , 1979; P r a k a s h , 1982). Hence the new analysis distinguishes the non-reactive degassing of such species as carbon dioxide or nitrogen from truly reactive phenomena in which the losses are intrinsic to the chemistry and microbiology of the subsurface e n v i r o n m e n t . T h e r e c h a r g e a n d d i s p e r s i v i t y a r e f u n c t i o n s o f t h e flow field, a n d accordingly may be used to characterize the related and important processes of aquifer reaeration and nutrient mixing with depth below the water table. These
0169-7722/89/$03.50
© 1989 Elsevier Science Publishers B.V.
94 latter phenomena are of interest to investigators of subsurface microbiology and in situ hazardous waste treatment. Our model builds upon the existing analytical modeling efforts of Gelhar and Wilson (1974) and Ostendorf et al. (1984), who derive the recharge hydraulics and source term behavior for the present study, as discussed below. The advectire-dispersion solutions of Ogata and Banks (1961) and Parker and Van Genuchten (1984) describe the open- and closed-system profiles respectively, although the transport is horizontal in the original contexts. Additional, less directly related, depth-integrated models may be cited as well: Lenau (1972) postulates a steady-state, conservative injection of pollution from a recharge well in a confined aquifer, while Wilson and Miller (1978) consider the transient response of a simply reactive contaminant from a similar source configuration. Chen (1987) derives a model including radial hydraulics and transport of a conservative species in a confined aquifer. Bear (1979) summarizes unsteady contaminant migration due to a series of one-dimensional reactive sources, and Prakash (1982) models steady-state reactive pollution due to instantaneous point, line, and volume sources. The method of characteristics adopted here follows one-dimensional analytical applications in the unsaturated and saturated zones by Wilson and Gelhar (1981) and Charbeneau (1981), respectively. Valocchi (1985, 1986) examines the validity of the local equilibrium assumption for simply adsorbing solutes in one-dimensional and radial flow fields on an analytical basis as well. The present simple analytical modeling approach complements existing more complex numerical contaminant transport codes. In this regard, we may distinguish two (Konikow and Bredehoeft, 1978; Borden and Bedient, 1986) and three (Huyakorn et al., 1986) -dimensional numerical solute transoirt models of single and coupled (Cederberg et al., 1985) constituents. These site-specific programs, with their attendant documentation requirements, are appropriate for the detailed analysis of fully measured plumes in complex geological settings. As such, they contrast with the present analytical effort, which is intended to provide a readily used, generic account of two-dimensional contamination of an idealized or sparsely measured flow field by a single species. HYDRAULICS We begin by considering the two-dimensional conservation of water mass in the unconfined aquifer, wherein a small (near)-vertical, average linear velocity component w is imposed by recharge upon a vertically uniform, (near)horizontal, average linear velocity u (Gelhar and Wilson, 1974): ~u ~x
--
+
0w ~z
-
0
(1)
with distance z below the water table and downgradient distance x from the source of pollution, as sketched in Fig. 1. Following Ostendorf et al. (1984), the unconfined aquifer flow field is perturbed by small changes due to a gently
95 landfill
x x x
x×x
It ~,w hs L _ _ ] -
~
source plane
----,
I
]
x,u h
Fig. 1. Definition sketch. sloping aquiclude, head loss, and recharge:
[ X(n
u = us 1 + ~ tanfl = tanfl'
s
-
tanfl
)1
nvus
kg
(2a) (2b)
with aquifer depth h, kinematic viscosity v, permeability k, gravitational acceleration g, recharge e, and porosity n. The s subscript denotes pollutant source conditions at the downgradient boundary of the landfill, while the slope tan fl relative to the water table includes the horizontally based slope tan fl'. The corresponding vertical average linear velocity follows from eqns. (1) and (2a). To leading order, we have: w -
n
- us tan
hs
(3)
Equation (3) suggests that w equals e / n at the water table and yields a velocity parallel to the aquielude where z equals ha, the first-order estimate of the actual aquifer thickness h. We are interested in frames of reference moving at speed u from the downgradient end of the landfill, taken as the (far field) source of pollution, as indicated by Fig. 1. The velocity and path xf of the frames are given by Ostendorf et al. (1984) as: dxf dt
--
(4a)
U
xf[
--tan )l
xf(~
z = -- 1 -U~ ~ss
rtUs
(4b)
96 =
t - ts
(4c)
with time r in the m o v i n g r e f e r e n c e frame which leaves the s o u r c e at time t~. A (near field) l i n e a r r e s e r v o i r analysis u n d e r the landfill has been derived by O s t e n d o r f et al. (1984) in o r d e r to specify s o u r c e plane c o n c e n t r a t i o n as a f u n c t i o n of ts and user population. Following the c o n c e p t of G e l h a r and Wilson (1974), this simple a p p r o a c h t r e a t s the complex n e a r field flow region as a fully mixed c h e m o s t a t r o u t i n g pollution input from the landfill to the o u t p u t plane at the d o w n g r a d i e n t end of the facility, w h e r e it serves as a v e r t i c a l l y uniform source of c o n t a m i n a t i o n . Thus, the frames c a r r y k n o w n pollution levels into the far field with presumed uniform initial v e r t i c a l profiles. OPEN- AND CLOSED-SYSTEM CONCENTRATION PROFILES The c o n s e r v a t i o n of n o n r e a c t i v e c o n t a m i n a n t mass c o n c e n t r a t i o n c can be c o n v e n i e n t l y studied in the m o v i n g r e f e r e n c e frame as a b a l a n c e of storage, v e r t i c a l advection, and v e r t i c a l dispersion: e~c ~ ~c ~2c ~-~ + -n--~z - D--0z2 = 0
(5a)
D
(5b)
= ~us
L o n g i t u d i n a l dispersion is neglected in this t r a n s p o r t e q u a t i o n on the premise of an appreciable l o n g i t u d i n a l velocity and a c o n t i n u o u s pollution source ( O s t e n d o r f et al., 1984); ~ a c c o r d i n g l y r e p r e s e n t s the v e r t i c a l dispersivity (Bear, 1979). The use of sin in place of w in the advective t r a n s p o r t term reflects o u r i n t e r e s t in c o n d i t i o n s n e a r the w a t e r table [see eq. (3)], w h e r e the bulk of t r a n s p o r t o c c u r s in the form of dilution and degassing t h r o u g h the r e c h a r g e lens. The g o v e r n i n g c o n s e r v a t i o n e q u a t i o n (5a) is solved subject to a presumedly uniform s o u r c e p l a n e c o n c e n t r a t i o n cs and a deep aquifer, i.e.: c =
cs
(T = 0)
(6a)
c =
cs
(z =
(6b)
~)
The l a t t e r c o n s t r a i n t suggests t h a t most of the t r a n s p o r t occurs n e a r the free surface and is u n a f f e c t e d by the bottom; an assumption in keeping with the use of E/n as the v e r t i c a l v e l o c i t y estimate. The r e m a i n i n g w a t e r table b o u n d a r y c o n d i t i o n differs for open and closed chemical systems: c =
ca
0c ~ - Dc~--~ + -nC = -nCa
(z = 0, open)
(7a)
(z = 0, closed)
(7b)
with c o n s t a n t a m b i e n t c o n c e n t r a t i o n c~. T h e open-system c o n d i t i o n (7a) rests on the premise of a rapidly diffusing u n s a t u r a t e d zone, as discussed for i n o r g a n i c c a r b o n in the Appendix.
97 Equations (5a), (6) and (7a) constitute the classical advective-diffusion problem with solution (Ogata and Banks, 1961): c~
ca
2
erfc L 2 ~ - - ~ J
+ exp ~-~ erfc L 2(Dr),. 2 j j
(8a)
(open) erfc y -
it,/2
exp ( - y ' 2 ) d y '
(8b)
Y
with complementary error function erfc y (Abramowitz and Stegun, 1972). Asymptotic behavior of this function may be used to advantage, whence: (z - ~ / n ) 2 7
c - Ca -- 1 -- i { rz__ ~,/n 1 (_~y,2exp ca - ca ~ erfcL2(Dr)l/2_j + 2
z + ~r/n
-4D-~T] ]1 z + ~ j
> 5(Dr) '/2
(9a) (9b)
The closed-system condition (7b) is predicated on a no-flux constraint at the water table, so that the advective and dispersive mechanisms are in balance at that location. A more complicated solution follows, as derived by Parker and Van Genuchten (1984):
c- Ca-
] + exp(~'Z )
cs -- ca
[Z +---~.1[~
~
( ~ ~1,,~ n \~)
[
erfc L 2(Dr)1/2 J
(z - ~/n) ~]
exp
]~-
+
1
~(z + 8"In)1 2On ] [~_ : ~/n l
-j - ~ erfc L 2(Dr) '/2 J
(10)
This too simplifies for a large argument: (z - ~r/n)'-']
C--Ca --
1+
c~ - ca
z + ~r/n
JlerfcrZ-
exp
z + ~:~/n > 5(Dr) 1/2
2
/n 1
(lea)
L2(-D~-~J (llb)
Figure 2 displays typical profiles for open and closed systems under dominant dispersion and recharge conditions for a pure ambient aquifer. The strong dispersion case features deeper mixing of open constituents and an appreciable free-surface closed-system concentration. Strong recharge on the other hand, establishes a pure overlying lens of water regardless of the open or closed nature of the pollutant. The area under the contaminant profiles is proportional to the depth averaged concentration: Fig. 2a indicates that open and closed systems exhibit values below their source levels. The open-system decrease is biggest since it includes dilution and degassing effects. The dilution reduction is predominant in Fig. 2b and is common to both systems.
98 0
~"--" ope n
-,~sed
0
"~, ~
\ open and closed
Z
hs
Z
0.5
hs 0.5
I
(°)
I.O o
c
o.5 / cs
,.o (b)
I
0
0,5
I.O
1.0
c~ cs
Fig. 2. Typical open and closed c o n t a m i n a n t system profiles for vs = 10-Sms -1, n = 0.3, and = 10Ss. S t r o n g dispersion (e = 1 0 - g m s -1, a = 10-~m) is sketched in Fig. 2a, while s t r o n g r e c h a r g e (~ = 10-Sins -~, ~ = 10-3m) is s h o w n in Fig. 2b. CALIBRATION AND T E S T I N G
The total inorganic carbon and chloride observations of Kimmel and Braids (1974, 1980), supplemented by raw data supplied by the USGS may be used to calibrate and test the open and closed models, respectively. Vertical and longitudinal profiles of contamination from a sanitary landfill on Babylon, Long Island provide ideal and independent data bases, since the hydraulics are well understood and the unconfined aquifer geometry is appropriately simple. Indeed, the existing depth-integrated model of Ostendorf et al. (1984) is calibrated with this plume, based on a source average linear velocity (us) of 3.37 × 10 6ms-l, recharge (e) of 3.25 × 10 9ms 1, aquifer depth (ha) of 22.5m at the source plane, permeability of 6.34 × 10-11m ~, viscosity of 1.1 × 10-6m2s -1, and a porosity of 0.27. Source plane concentrations are computed from monitoring well data located at the downgradient boundary of the landfill, following the near field procedure proposed by Ostendorf et al. (1984). The source plane wells are used to calibrate a per capita contaminant generation rate Scl equals 1.40 × 10 _8 kg c a p - i s -1 for chloride in the earlier analysis, and this value is adopted here to estimate the ca values in closed system testing. A total inorganic carbon value of Sc equals 7.90 × 10-gkgcap-ls -1 generates the open-system source term concentrations in the present study, based on 1973 observations in the source plane monitoring wells. The new Sc value is needed since Ostendorf et al. (1984) studied bicarbonate and not total inorganic carbon in their prior work. We compare data and theory using statistics of the error 5 defined by: 5 =
c (measured) - c (predicted) c (measured)
(12a)
= _1.Z5 ]
(12b)
= ~ . Z&2 - ~ t I/2
(12c)
99
T h e m e a n e r r o r ~ and s t a n d a r d d e v i a t i o n a are c o m p u t e d in a c c o r d a n c e with B e n j a m i n and Cornell (1970). The sign of 3 indicates model over or underprediction and is a c c o r d i n g l y useful in identifying systematic model errors. The e r r o r s t a n d a r d d e v i a t i o n is based on the absolute v a l u e of individual 3's and conseq u e n t l y m e a s u r e s the m a g n i t u d e of the error. In this regard, a b o u t 2/3 of our predictions lie w i t h i n a of t h e i r m e a s u r e d values for a zero m e a n error. T h e (far field) l o n g i t u d i n a l total i n o r g a n i c c a r b o n profile d a t a observed in 1973 by Kimmel and Braids (1980) are used to c a l i b r a t e the vertical dispersivity in eqns. (5b) and (8a) with the results summarized in Table 1 and Fig. 3a. The TABLE 1 L o n g i t u d i n a l profile of total inorganic carbon: model calibration x (m)
~ (10 s s)
cs (kg m - 3)
z (m)
c (meas) (kg m 3)
c (pred) (kg m - 3)
(%)
360
1.05
0.136
630 900
1.80 2.54
0.117 0.091
920
2.59
0.089
1550 1890 2180
3.89 5.02 5.69
0.057 0.045 0.038
2230 2260 2810
5.80 5.87 7.05
0.037 0.036 0.020
1.75 9.05 25.11 23.16 3.08 8.97 18.45 11.74 21.50 12.17 12.69 8.75 19.33 17.66 21.38 11.89 18.26
0.017 0.087 0.141 0.091 0.007 0.136 0.137 0.035 0.060 0.055 0.067 0.018 0.037 0.010 0.040 0.010 0.010
0.022 0.100 0.136 0.114 0.015 0.044 0.077 0.055 0.081 0.029 0.020 0.011 0.024 0.021 0.024 0.007 0.010
- 32 - 15 4 24 - 115 68 44 56 - 36 47 69 42 36 99 40 32 - 6
Notes: 1973 data (Kimmel and Braids, 1980), supplemented by USGS r a w data. Source conditions from O s t e n d o r f et al. (1984). C o n c e n t r a t i o n s expressed as carbon. Calibrated dispersivity ~ 0.073 m, zeros m e a n error. Calibration e r r o r s t a n d a r d deviation a = 54%. 0.15
el
I
0.10 "u'x"~. ~
c (kg m-3)
I I0
"a = 0.073 m
Z (rn)
o
0.05__
0
I
©15-25m o5-15m eO-5m
(0) 0
II
20
(b) I xaXt
50 500
1000 x
1500 (m)
2000
2500
3000
0
l
0.025 0.050 0.0"75 ¢(kg m -3)
0.0[0
Fig. 3. Total c a r b o n calibration (a) and test (b) data (circles) and predictions (lines). 1973 observations of Kimmel and Braids (1974, 1980), supplemented by raw data supplied by USGS.
100
more distant wells have larger travel times (~), corresponding to frames leaving the landfill at earlier times with smaller user populations hence lower source concentrations. The vertical location refers to distance below the average water table elevation, and the open system is chosen for calibration purposes since it possesses a steeper gradient than its closed count erpart and is accordingly more sensitive to ~. The calibration is a Fibonacci search (Beveridge and Schechter, 1970) for the dispersivity value zeroing the mean error. The operation yields a standard deviation of 54% associated with the optimal dispersivity: = 0.073 m
(13)
This dispersivity value is close to the first-order decay based value of 0.10 m put forth by Ostendorf et al. (1984) on an ad hoc basis in their depth-integrated analysis of bicarbonate data at the site. The newer version reflects a more rigorous estimate of the transport mechanisms: the "decay" is in reality a mass transport process constrained by vertical mixing and recharge. Table 1 suggests no systematic dependence of the error upon horizontal or vertical distance, so that the relatively high standard deviation may in part be attributed to sampling errors. In any event, additional independent testing of the calibrated model is certainly warranted. The vertical inorganic carbon, longitudinal chloride, and vertical chloride observations of Kimmel and Braids (1974, 1980) provide these necessary tests, as summarized in Tables 2-4 and Figs. 3b and 4. The accuracy of the vertical inorganic carbon test is encouraging, as evidenced by the 28% mean and 21% standard deviation error statistics. The longitudinal chloride data provide an additional measure of support, particularly if one neglects two data points immediately beneath the water table. Most of the mean error ( - 4 6 % ) and standard deviation (126%) for this test is created by the two near-surface data points 360m and 900m downgradient: we compute quite accurate values of - 2°./o and 25% for the statistics without these samples. The vertical chloride TABLE 2 Vertical profile of total inorganic carbon: model test z (m)
c (measured) ( k g m ~)
c (predicted) ( k g m -3)
3 (%)
5.58 12.01 14.27 19.30 24.17
0.0152 0.0382 0.0405 0.0921 0.0420
0.0113 0.0249 0.0294 0.0378 0.0442
25 34 27 59 - 5
Notes: USGS r a w data, except 24.17 m value (Kimmel and Braids, 1974). Source condition from O s t e n d o r f et al. (1984). C o n c e n t r a t i o n s expressed as carbon. Profile location x = 1570m, r = 4.26 × 108s. Test statistics ~ = 28%, a = 21%.
101 TABLE 3 L o n g i t u d i n a l profile of chloride: model t e s t x (m)
~ (10Ss)
c~ ( k g m 3)
z (m)
c (meas) ( k g m -3)
c (pred) ( k g m 3)
(%)
360
1.05
0.258
630 900
1.80 2.54
0.224 0.177
920
2.59
0.174
1550 1890 2180
3.89 5.02 5.69
0.117 0.096 0.087
2230 2260 2810
5.80 5.87 7.05
0.080 0.079 0.051
1.75 9.05 25.11 23.16 3.08 8.97 18.45 11.74 21.50 12.17 12.69 8.75 19.33 17.66 21.38 11.89 18.26
0.054 0.190 0.270 0.180 0.024 0.200 0.360 0.130 0.180 0.109 0.150 0.060 0.057 0.057 0.067 0.042 0.045
0.213 0,247 0.258 0.223 0.134 0.156 0.172 0.159 0.171 0.102 0.081 0.067 0.078 0.071 0.072 0.042 0.045
- 294* - 30 5 - 24 - 459* 22 52 - 22 5 6 46 - 12 - 38 - 25 - 8 0 0
Notes: 1973 d a t a ( K i m m e l a n d Braids, 1980), s u p p l e m e n t e d by U S G S r a w data. S o u r c e c o n d i t i o n s from O s t e n d o r f et al. (1984). A m b i e n t c o n c e n t r a t i o n of 0.020 k g / m 3 included. T e s t s t a t i s t i c s w i t h n e a r - s u r f a c e p o i n t s ~ = - 46%, ~ = 126%. * T e s t s t a t i s t i c s w i t h o u t n e a r - s u r f a c e p o i n t s ~ = - 2%, ~ = 25%.
TABLE 4 V e r t i c a l c h l o r i d e profile: model t e s t z (m)
c (measured) ( k g m -3)
c (predicted) ( k g m -3)
(%)
5.58 12.01 14.27 19.30 24.17
0.088 0.140 0.170 0.200 0.170
0.082 0.094 0.097 0.102 0.106
6 33 43 49 39
Notes: U S G S r a w data, e x c e p t 24.17m d a t a ( K i m m e l a n d Braids, 1974). S o u r c e c o n d i t i o n s from O s t e n d o r f et al. (1984). Profile l o c a t i o n x = 1570m, ~ = 4.26 × 108s. T e s t s t a t i s t i c s ~ ~ 34%, ~ = 15%. test summarized the
34%
ambient
mean chloride
the plume.
in Table and
15%
4 and Fig. 4b does not fare as well however, standard
concentration
The closed-system
deviation.
of 0.020 kgm errors
do suggest
Both
chloride
-3, as indicated a tendency
tests
based
on
reflect
an
by wells
to overpredict
outside near
102 0.40
I
01
I
I
,,5!25m o5-15m •O-5m
0.30 c (kg m-3)
o % .
0.20
0
I
ct= 0.073 m
%.
- -
I0
Z (m)
2O
0. I0
(o)
I 500
0
"i I000
1
I
1500
2000
I
(b)
30
2500
0
3000
I
0.05
I
O.lO 0.15 c (kg m-3)
0.20
X (m]
Fig. 4. Longitudinal (a) and vertical (b) chloride test data (circles) and predictions (lines). 1973 observations of Kimmel and Braids (1974, 1980), supplemented by raw data supplied by USGS. the free surface and u n d e r p r e d i c t at d e e p e r locations in the profile, perhaps in response to a n o n u n i f o r m s o u r c e b e h a v i o r not included in the present analysis. T a k e n as a whole then, the c a l i b r a t i o n and test results provide a modest level of s u p p o r t for this analysis, a l t h o u g h a d e p t h - v a r y i n g s o u r c e term c h a r a c t e r i z a tion could well improve the model a c c u r a c y for the closed-system constituents. RESULTS T h e c a l i b r a t e d and tested open-system model m a y be c o m p a r e d with its one-dimensional c o u n t e r p a r t to e l u c i d a t e the physical basis of the ad hoc decay c o n s t a n t A used in c o n v e n t i o n a l models. T h e d e p t h - a v e r a g e d c o n c e n t r a t i o n follows upon i n t e g r a t i o n of eqn. (8a), as derived in the Appendix: ca cs - c.
1
ierfc y
~
=
h,
ierfc
erfc(y')dy"
':2]
-
erf
k
]
(14a)
(14b)
Y
erfy
=
1 - ierfcy
(14c)
with e r r o r f u n c t i o n i n t e g r a l ierfc y defined by A b r a m o w i t z and S t e g u n (1972). E q u a t i o n (14) m a y be used with d e p t h - i n t e g r a t e d d a t a to c a l i b r a t e a vertical dispersivity v a l u e in the absence of m o r e definitive v e r t i c a l open-system conc e n t r a t i o n profiles. T h e " t r u e " one-dimensional model r e p r e s e n t e d by eqn. (14) c o n t r a s t s with the c o n v e n t i o n a l b a l a n c e of a d v e c t i o n and first-order d e c a y (Bear, 1979): - ca =
(cs - c , ) e x p ( - / l ~ )
(15)
in t h a t r e c h a r g e and v e r t i c a l dispersivity p a r a m e t e r s a p p e a r explicitly in the f o r m e r e x p r e s s i o n as a c o n s e q u e n c e of the p r i o r d e p t h - v a r y i n g analysis. The ad
103
hoc decay constant can be related directly to the responsible vertical transport properties for small ~ values, whence: --
Ca
Cs - -
Ca
(1 - ~tr)
- C, cs - Ca
1 -
2 I__~l I/2 ~
(ad hoc decay)
(16a)
(true mixing)
(16b)
Thus dispersion induces pseudoreactive decay of a degassing constituent, with an ad hoc constant approximately given by: -
)+ ~
(17a)
-
h+ \ ~ X c /
xc = u+~
(17b)
with the plume length scale Xc in the longitudinal direction specified by advection in the unconfined aquifer. We may also estimate the flux J of open system contaminant through the water table once the concentration profile is specified, since J is the combined result of vertical recharge advection and dispersion (Bear, 1979): J
=
nw(c-
c.) -
Oc
nD-~z
(z = 0)
(18)
In view of eqns. (7a) and (8), we find:
g = -(c+-
{i/D'~l/2 [- (~.c/n)2_l
ca) n[~)__
exp|_
4D+ J -
~
-/"
~.1/2
x)
~erfc[~)~
(19)
The contrasting roles of the vertical transport mechanisms may be noted for limiting cases of strong recharge and strong dispersion J
~
-n(c~
J
~ 0
(O~']2
- c . ) \-~-~/
=
~)
(20a)
(s =
~)
(20b)
(D
Dispersion mixes the carbon towards the water table at a rate that decreases with travel time due to a flattening vertical gradient. The vertical recharge flow flushes pollution back into the plume by downward advection, and so reduces J. CONCLUSIONS
We derive simple two-dimensional mathematical models describing the unsteady transport of conservative contaminants through an unconfined aquifer with a gently sloping aquiclude subject to advection, recharge, and vertical dispersion. Open (dispersive) and closed (nondispersive) chemical
104 systems are considered by imposition of appropriate boundary conditions at the water table, and the pseudoreactive depth-integrated behavior exhibited by the volatile open-system constituents is explained in terms of fundamental vertical mixing parameters. Total inorganic carbon and chloride data at the Babylon, Long Island landfill are used to calibrate and test the open- and closed-system models with a modest level of accuracy. The observed longitudinal distribution of total inorganic carbon is used for calibration, an exercise yielding a vertical dispersivity estimate of 0.073m with a standard deviation of 54%. Vertical profiles of inorganic carbon and chloride as well as a longitudinal chloride data set are used to test the theory, an exercise resulting in means and standard deviations ranging from - 2 to 34% in magnitude. The statistics do not include two chloride data points immediately below the water table, which are substantially overpredicted by the model. Future research may proceed along various related fronts. The present analysis of groundwater contamination rests on a sequential understanding of the following phenomena, both in terms of data and theory: - - groundwater hydraulics - - pollution source history depth-integrated contamination, conservative contaminants depth-varying contamination, conservative contaminants The first extension of the modeling effort should be the inclusion of a depthvarying source condition in response to the systematic errors noted in chloride testing. Such a modification should preserve the simplicity inherent in an analytical model approach while maintaining the accuracy of the open-system profile tests in the present study. Additional extensions of the modeling effort might include the depth-varying transport of truly reactive pollutants, first for a single species and then for coupled contaminants. One is certainly drawn to dissolved oxygen as a reactive constituent of interest, and the speciation of inorganic carbon at Babylon seems a logical starting point for the coupled work. The present model should be modified to account for the influence of a nondispersive bottom boundary to expand its applicability, with the possible inclusion of an appropriately simple allowance for the vertical variation of the vertical velocity w. Additional testing of the present theory is certainly called for, at other sites and for other contaminants. In the latter regard, total nitrogen in a sewage effluent plume seems a likely candidate for open system testing, due to degassing associated with the denitrification process. -
-
-
-
ACKNOWLEDGMENTS Partial funding for this research was provided by the Massachusetts Division of Water Pollution Control under Grant Number 86-30765; the Authors acknowledge and appreciate their support. We also thank the Water Resources Division of the United States Geological Survey for providing the raw data base for the Babylon landfill.
105 REFERENCES Abramowitz, M. and Stegun, I.A., 1972. Handbook of Mathematical Functions. National Bureau of Standards, Washington, DC, 1046 pp. Bear, J., 1979. Hydraulics of Groundwater. Wiley, New York, NY, 567 pp. Benjamin, J.R. and Cornell, C.A., 1970. Probability, Statistics, and Decision for Civil Engineers. McGraw-Hill, New York, NY, 684 pp. Beveridge, G.S.G. and Schechter, R.S., 1970. Optimization: Theory and Practice. McGraw-Hill, New York, NY, 773 pp. Borden, R.C. and Bedient, P.B., 1986. Transport of dissolved hydrocarbons influenced by oxygen limited biodegradation 1. Theoretical development. Water Resour. Res., 22: 1973-1982. Cederberg, G.A., Street, R.L. and Leckie, J.O., 1985. A groundwater mass transport and equilibrium chemistry model for multicomponent systems. Water Resour. Res., 21:1095 1104. Charbeneau, R.L., 1981. Groundwater contaminant transport with adsorption and ion exchange chemistry: method of characteristics for the case without dispersion. Water Resour. Res., 17: 705 713. Chen, C.S., 1987. Analytical solutions for radial dispersion with Cauchy boundary at injection well. Water Resour. Res., 23:1217 1224. Gelhar, L.W. and Wilson, J.L., 1974. Groundwater quality modeling. Ground Water, 12: 399~408. Hillel, D., 1982. Introduction to Soil Physics. Academic Press, New York, NY, 364pp. Huyakorn, P.S., Jones, B.G. and Andersen, P.F., 1986. Finite element algorithms for simulating three-dimensional groundwater flow and solute transport in multilayer systems. Water Resour. Res., 22:361 374. Kimmel, G.E. and Braids, O.C., 1974. Leachate plumes in a highly permeable aquifer. Ground Water, 12: 388-393. Kimmel, G.E. and Braids, O.C., 1980. Leachate plumes in groundwater from Babylon and Islip landfills, Long Island, New York. U.S. Geol. Surv., Prof. Pap. 1085, 38 pp. Konikow, L.F. and Bredehoeft, J.D., 1978. Computer model of two-dimensional solute transport and dispersion in groundwater. Techn. Water Resour. Invest. USGS, USGS, Reston, VA. Lenau, C,W., 1972. Dispersion from recharge well. J. Eng. Mech. Div., ASCE, 98: 331-334. Ogata, A. and Banks, R.B., 1961. A solution of the differential equation of longitudinal dispersion in porous media. U.S. Geol. Surv., Prof. Pap. 411-A, 7 pp. Ostendorf, D.W., Noss, R.R. and Lederer, D.O., 1984. Landfill leachate migration through shallow unconfined aquifers. Water Resour. Res., 20: 291-296. Parker, J.C. and Van Genuchten, M.T., 1984. Flux averaged and volume averaged concentrations in continuum approaches to solute transport. Water Resour. Res., 20: 866-872. Prakash, A., 1982. Groundwater contamination due to vanishing and finite size continuous sources. J. Hydraul. Div., ASCE, 108: 572-590. Stumm, W. and Morgan, J.J., 1981. Aquatic Chemistry. Wiley, New York, NY, 780 pp. Valocchi, A.J., 1985. Validity of the local equilibrium assumption for modeling sorbing solute transport through homogeneous soils. Water Resour. Res., 21: 808-820. Valocchi, A.J., 1986. Effect of radial flow on deviations from local equilibrium during sorbing solute transport through homogeneous soils. Water Resour. Res., 22: 1693-1701. Wilson, J.L. and Gelhar, L.W., 1981. Analysis of longitudinal dispersion in unsaturated flow 1. The analytical method. Water Resour. Res., 17:122 130. Wilson, J.L. and Miller, P.J., 1978. Two-dimensional plume in uniform groundwater flow. J. Hydraul. Div., ASCE, 104: 503-514. APPENDICES
Inorganic carbon boundary condition T h e i n o r g a n i c c a r b o n c o n c e n t r a t i o n c o n s i s t s o f b i c a r b o n a t e cl a n d d i s s o l v e d c a r b o n d i o x i d e c2 c o n t r i b u t i o n s ( S t u m m a n d M o r g a n , 1981) i n t h e p H r a n g e
106
(5-7) common to sanitary landfill plumes like t h a t at Babylon, Long Island: c = cl + c2
(21a)
cl = c2 ~
(21b)
K~ = 10-63mol L -1
(21c)
The equilibrium equation (21b) includes the hydrogen ion activity [H +] and equilibrium constant K1 expressed in moles/liter. The latter parameter is temperature and concentration dependent to a modest degree, but our intent here is an order of magnitude discussion of mass fluxes and a representative value is sufficient. The inorganic carbon and its species are expressed in terms of k g m -3 carbon, for consistency with the main analysis. The carbon dioxide species governs the total inorganic carbon degassing rate, since it is in gaseous equilibrium with the unsaturated zone air monolayer at the water table. Henry's law is appropriate (Stumm and Morgan, 1981) at this location, and may be expressed as:
c2 = KHRTp
(z = 0)
(22a)
R = 0.0821 atm L m o l - ' K 1
(22b)
KH = 0.045 mol L - ' a t m - '
(22c)
with Henry's constant KH, universal gas constant R, temperature T, and gaseous carbon dioxide density p, in k g m -3 carbon. For typical groundwater temperatures of 290 K, eqn. (22) suggests that the dissolved concentration and gaseous density are of equal magnitude on either side of the free surface: c2 ~
p
(z = 0)
(23)
The carbon dioxide (hence total inorganic carbon) flux J through the unsaturated zone is a gaseous dispersion process which may be estimated in accordance with (Hillel, 1982): J
~
-D~ P
(24)
with unsaturated zone thickness b and gaseous dispersivity Da of order 10- 6 m2s-, in magnitude. This crude estimate may be compared with its liquid side counterpart: J
~
c
-D-h~
(25)
in order to assess the magnitude of c at the water table. Collecting eqns. (21)-(25), we deduce: c(z = 0) c
D(
D~ 1 + ~
K1)
(26a)
107 C(Z = 0) ~
C
(pH
< 7.3)
(26b)
where the s a t u r a t e d and u n s a t u r a t e d zone thicknesses are the same order of m a g n i t u d e in size. P h y s i c a l l y speaking, the u n s a t u r a t e d zone efficiently carries off w h a t e v e r total carbon the s a t u r a t e d zone can deliver due to the smallness (1/10) of the dispersivity ratio; the simple b o u n d a r y condition (7a) is the result.
Open system The depth-integrated open-system c o n c e n t r a t i o n is defined as: hs
1 f cdz
-
(27)
hs
so that, recalling eqn. (8a): = c~ I =
i
cs - ca 2hs (I + II)
(28a)
L[zb/--n~--~/dz - ~ln]
(28b)
erfc
0 II --
i (sn)~ exp
L[z+~z/nl-2~D~-~ J
erfc . - d z
(28c)
0
The use of the infinite upper limit of i n t e g r a t i o n in eqn. (28) reflects the assumption of an infinite bottom b o u n d a r y condition in eqn. (6b). The first integral I is straightforward, while the second II requires an i n t e g r a t i o n by parts: II
nD~erfc
-
II -
+
nD [ erfc ( n i p2=~''cli2"~ ~, ]i ~ . -
(D~ 'i2
\~--~]
erfc
i[ exp
-_ aUn)~.~ 4D~
\2nD'121#_] (e'rli2"~]
j dz
(29a)
(29b)
This last relation simplifies with the aid of an identity: erfc(-y)
= 2 - erfcy
giving rise to eqn. (14a).
(30)