- Email: [email protected]
Geostatistical schemes for groundwater sampling
Journal of Hydrology, 103 (1988) 85 102 Elsevier Science Publishers B.V., Amsterdam
85 Printed in The Netherlands
[4] G E O S T A T I S T I C A L S C H E M E S FOR G R O U N D W A T E R S A M P L I N G
SHAHROKH ROUHANI and TIMOTHY J. HALL
School of Civil Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332 (U.S.A.) (Received March 12, 1987; revised and accepted January 8, 1988)
ABSTRACT Rouhani, S. and Hall, T.J., 1988. Geostatistical schemes for groundwater sampling. J. Hydrol., 103:85 102. Geostatistical techniques offer efficient tools for design of groundwater sampling networks. They include procedures for the selection of the best sequences of sampling points, such as: variance reduction analysis, median ranking, and risk ranking. Variance reduction analysis considers primarily the accuracy of the estimated field, while median ranking is based only on the magnitude of the estimated values. Risk ranking is a compromise between these procedures that appears to yield a more balanced guideline for cases when planners desire to acquire maximum information, while monitoring areas where the variable of interest exhibits critical values. These procedures are used for the design of a regional shallow groundwater quality monitoring network in the Dougherty Plain, located in southwest Georgia. The shallow aquifer of concern is the main recharge route to a semiconfined aquifer which is the primary source of water in this region. The desired monitoring network acts as an early warning system for groundwater po]lution in deeper layers. Leakance data are utilized to identify the primary sampling locations. Statistical analyses indicate that leakance has a log-normal distribution with a constant drift and a linear spatial covariance. The results of our risk rankings demonstrate that the southern tip of the Dougherty Plain and its upper central zone should be the prime targets of our monitoring activities.
INTRODUCTION
An increasing number of researchers are advocating the use of geostatistical techniques in the design of groundwater sampling networks (Rouhani and Fiering, 1986). For this purpose the sampled variables should be viewed as stochastic processes. In particular, spatial variables such as the transmissivity, storativity, and steady-state piezometric heads may be considered as random fields, with the following structure: Z(X,)
=
M(X~) + R ( ~ )
(1)
where: Z(X/) = the actual value of the variable of interest at X/; M(X/) = a slowly varying deterministic function known as the drift which is equal to the expected value of Z at X~; and R(X~) = a spatially fluctuating random function with zero expectation and a spatial covariance function. As Rouhani (1985) states, the data management of geohydrologic random
0022-1694/88/$03.50
© 1988 Elsevier Science Publishers B.V.
86 fields depends upon the objectives of the planning approach. Often very little data are available. Furthermore, the measured values may be clustered together and therefore not provide information about the whole field. In such situations, planners may wish to design a data collection scheme in order to better define their variables of interest. The questions that immediately arise are: How should we design such a regional sampling network? More specifically, what criteria should be utilized as the basis of our network design? Where are the best locations for sampling sites? Do the selected sites remain unchanged if we change our sampling criterion? Geostatistical techniques such as kriging provide tools to answer these questions. Commonly geostatistical sampling procedures are based on the maximization of incremental information subject to budget constraints. For instance, Rouhani (1985) proposes variance reduction analysis to select those sequences of n points so chosen from m potential sites to maximize reduction in the total variance of estimates. In these procedures the criteria of site selection are such that they give more priority to points with high estimation variances, regardless of their estimated magnitudes. These criteria are suitable for cases where the estimated value of the variable of concern is not of primary importance. In many cases, however, planners are not only interested in gaining as much information as possible, but also to monitor areas where the sampled variable exhibits critical values. This requires some modifications in the selection criterion to include both the accuracy and the magnitude of the estimated values. To accomplish this task we propose to utilize the risk value as an alternative selection criterion, as described in the following sections. GEOSTATIST1CALANALYSIS Geostatistical techniques provide means for the statistical analysis of spatially distributed random variables. Matheron (Journel and Huijbregts, 1978) states that: ~Geostatistics is the application of the formalism of random functions to the reconnaissance and estimation of natural phenomena." At first glance, the measured values of many geohydrological variables appear to be random and unpredictable. However, values taken at neighboring points may be related by statistical functions. Geostatistics enables us to estimate these statistical relations and utilize them to estimate the value of the variable at an arbitrary point. Kriging is one of the important geostatistical techniques which has been applied to many geohydrologic problems for estimation and sampling purposes. This algorithm is basically a statistical interpolation technique. In point kriging one estimates the value of a random field at any arbitrary point. X~), from values measured at other points, in a linear form of: N
Z(Xo) = ~ )~ioZ(X~) i
1
(2)
87
where: Z(X0) = k r i g i n g e s t i m a t e at 3/o; Z(Xi) - m e a s u r e d v a l u e at Xi, i - 1 ..... N; 2li0 - k r i g i n g w e i g h t for Z(X~) to e s t i m a t e Z(X0); a n d N - n u m b e r of measu r e m e n t points. T h e 2i0 are defined by two criteria: (1) unbiasedness; E[2(X0) Z(X~,)] - 0, w h e r e Z(Xo) is the a c t u a l v a l u e of the field at X0; and (2) m i n i m u m s q u a r e d e r r o r of estimation: this r e q u i r e s E[Z(X0) - Z(X0)] 2 to be m i n i m u m . The m i n i m u m s q u a r e d e r r o r of e s t i m a t i o n is also a m e a s u r e for the a c c u r a c y of o u r e s t i m a t e s which is k n o w n as e s t i m a t i o n (or kriging) v a r i a n c e . K r i g i n g a s s u m e s t h a t the v a r i a b l e of i n t e r e s t is a r a n d o m field w i t h a s t r u c t u r e s i m i l a r to eqn. (1). In this w o r k w i t h o u t a n y loss of g e n e r a l i t y , it is a s s u m e d t h a t the e x p e c t e d v a l u e of Z(X0), the so-called drift, is a p o l y n o m i a l of kth order: M(Xo)
l(k) ~ bJp(Xo)
, p
(3)
1
where: b t , - fixed u n k n o w n coefficients; f p ( X ) = the p t h m o n o m i a l ; and l(k) = the n u m b e r of m o n o m i a l s in the a b o v e kth o r d e r polynomial. In a t w o - d i m e n s i o n a l s p a c e with C a r t e s i a n c o o r d i n a t e s (x,y), a s e c o n d - o r d e r p o l y n o m i a l (k - 2, and 1 - 6) h a s the following m o n o m i a l s : 1, x, y, xy, x 2, and y2.
T h e r a n d o m c o m p o n e n t in eqn. (1), R, is a s s u m e d to h a v e a zero m e a n and a c o v a r i a n c e , defined as K[R(Xi),R(X~)] = E [ R ( X i ) R ( X j ) ]. It is f u r t h e r a s s u m e d t h a t the c o v a r i a n c e is an isotropic f u n c t i o n of the d i s t a n c e b e t w e e n a n y two points: K 0 = K([X/ - Xj[). T h e e s t i m a t i o n of the c o v a r i a n c e is a c c o m p l i s h e d u s i n g the s t r u c t u r a l a n a l y s i s p r o p o s e d by Delfiner (1975). The s u g g e s t e d p o l y n o m i a l c o v a r i a n c e f u n c t i o n h a s the following form: k
K(h)
Cb(h) +
~ p
av.
,h 2p'1
(4)
1
where: K - the c o v a r i a n c e function; h = the length of v e c t o r d i s t a n c e b e t w e e n a n y two points; ap - the p t h coefficient; C = the point v a r i a n c e ; = the D i r a c d e l t a function; and k = the o r d e r of the p o l y n o m i a l drift. T h e e s t i m a t i o n v a r i a n c e c a n t h e n be c a l c u l a t e d as:
Po i
N
N
2
2 0j
:
)'io~tjoK(IXi - XJt )
(5)
0
w h e r e Vo is the e s t i m a t i o n v a r i a n c e at X0, Zoo - - 1, and the r e s t of definitions are a l r e a d y given. F o r f u r t h e r i n f o r m a t i o n is r e f e r r e d to R o u h a n i (1983 and 1985). V a r i a n c e reduction a n a l y s i s T h e e s t i m a t i o n v a r i a n c e in eqn. (5) is a m e a s u r e for the a c c u r a c y of the e s t i m a t e d v a l u e at X0. Utilizing this function, v a r i a n c e r e d u c t i o n a n a l y s i s provides a s c h e m e based on the m a x i m i z a t i o n of i n c r e m e n t a l i n f o r m a t i o n
88 subject to b u d g e t c o n s t r a i n t s (Rouhani, 1985). This p r o c e d u r e first m e a s u r e s the m a g n i t u d e of i n f o r m a t i o n g a i n at a n y e s t i m a t e d point, X0, due to a m e a s u r e m e n t at X * (the a r b i t r a r y l o c a t i o n of a p o t e n t i a l s a m p l i n g site), in t e r m s of r e d u c t i o n s in V0 a f t e r to the new sampling, as: N
VR,~.
ltk)
[Ko. = ~ i
2i, K,o1
~ p
l,p.fp(Xo)['/[V.(N)]
(6)
1
where: VR.+ = the v a r i a n c e r e d u c t i o n (i.e., a c c u r a c y or i n f o r m a t i o n gain) at X¢~ due to a m e a s u r e m e n t at X,; Ki~ the c o v a r i a n c e f u n c t i o n b e t w e e n X, and Xj; 2,+ = the o p t i m a l w e i g h t of Z(Xi) in e s t i m a t i o n of Z(X.) p r i o r to the new m e a s u r e m e n t ; /tp, = the L a g r a n g e m u l t i p l i e r for the p t h m o n o m i a l filtering c o n s t r a i n t in the k r i g i n g s y s t e m for the e s t i m a t i o n of Z(X.) p r i o r to the new m e a s u r e m e n t ; N = the n u m b e r of existing d a t a points prior to the new meas u r e m e n t ; and V,(N) the k r i g i n g v a r i a n c e at X, p r i o r to the new m e a s u r e ment. It t h e n e x p a n d s the a b o v e m e a s u r e to c o v e r the e n t i r e field as follows:
TVR.
~ VRj,
(7)
i
where: TVR, = the t o t a l v a r i a n c e r e d u c t i o n due to a m e a s u r e m e n t at X,; and j - the set of e s t i m a t e d points. TVR, r e p r e s e n t s the a m o u n t of r e g i o n a l inf o r m a t i o n g a i n due to a new m e a s u r e m e n t at X,. At e a c h r o u n d of sampling, the site a m o n g p o t e n t i a l s a m p l i n g l o c a t i o n s with m a x i m a l TVR, is selected as the n e x t m e a s u r e m e n t location. This yields a s e q u e n c e of n points a m o n g m sites for f u r t h e r sampling.
Risk ranking T h e a b o v e a p p r o a c h gives m o r e p r i o r i t y to points with h i g h e s t i m a t i o n v a r i a n c e s , r e g a r d l e s s of t h e i r e s t i m a t e d m a g n i t u d e . S u c h a c r i t e r i o n is not s u i t a b l e for a case w h e n p l a n n e r s desire to g a i n m a x i m u m i n f o r m a t i o n , while m o n i t o r i n g a r e a s w h e r e the v a r i a b l e of c o n c e r n exhibits critical values. F o r i n s t a n c e , in o u r case study, the desired n e t w o r k is aimed to m o n i t o r the q u a l i t y of r e c h a r g e d w a t e r at critical l e a k y points. So, we are not o n l y i n t e r e s t e d in g a i n i n g as m u c h i n f o r m a t i o n as possible, b u t also to m o n i t o r a r e a s with p o t e n t i a l l y h i g h levels of r e c h a r g e . This implies t h a t the s e l e c t i o n c r i t e r i o n should be modified to include b o t h the a c c u r a c y a n d the m a g n i t u d e of o u r e s t i m a t e d values. T h i s o b j e c t i v e c a n be a c c o m p l i s h e d by utilizing the risk v a l u e as an altern a t i v e selection criterion. T h e risk v a l u e is defined as the v a l u e of the v a r i a b l e of i n t e r e s t w h o s e p r o b a b i l i t y of e x c e e d e n c e is ~ percent. F o r a l o g - n o r m a l l y d i s t r i b u t e d v a r i a b l e the risk v a l u e can be w r i t t e n as: Z~(X) =
exp[E{ln[Z(X)]} + z, (Var{ln[Z(X)l}) 1'~]
(8)
where: Z~(X) = the risk v a l u e at X w i t h the p r o b a b i l i t y of e x c e e d e n c e of a
89
.
ANTA
DOUGHERTY
Fig. 1. Geological profile of Georgia, U.S.A. percent; Z(X) = r a n d o m variable of interest at X; Var( ) = the variance; and z~ - the standardized n o r m a l l y distributed r a n d o m variable with a probability of e x c e e d e n c e of ~-percent. In o u r case study, at each r o u n d of kriging, the point with maximal risk value is identified as the n e x t best sampling site, which will yield an ordered list of points for site selections. As we r e d u c e the level of o u r risk, ~, more weight will be given to the variance. It can be inferred t h a t v a r i a n c e r e d u c t i o n analysis is an e x t r e m e case of risk r a n k i n g when v a r i a n c e has the complete weight. On the o t h e r extreme, when ~ = 50%, we h a v e the m e d i a n r a n k i n g which totally ignores the a c c u r a c y of the estimates. APPLICATION For o u r case study we have selected the g r o u n d w a t e r q u a l i t y m o n i t o r i n g in the shallow aquifer of the D o u g h e r t y Plain in s o u t h w e s t e r n Georgia, as shown in Fig. 1. This a r e a is a rapidly growing a g r i c u l t u r a l region, which is u n d e r l a i n by a succession of sand, clay and c a r b o n a t e rocks to a depth of more t h a n 5000 ft, forming one of the most p r o d u c t i v e m u l t i l a y e r aquifer systems in the United States. The a b u n d a n t supply of good quality g r o u n d w a t e r , along with a mild climate, and a flat to gently rolling t e r r a i n has induced significant g r o w t h in the a g r i c u l t u r a l p r o d u c t i o n of this region. This growth is accom-
90 panied by an increased use of fertilizers and pesticides. Recent data indicates that the average fertilizer use in the study area is approximately 105 t acre In 1977 this level was even higher at about 150 t acre l (Georgia Department of Agriculture, 1985). Some of the chemical components of these fertilizers and pesticides are toxic to human, long-lasting, and tend to accumulate in the hydrogeological system. In a recent study by the Georgia Department of Natural Resources, pesticides were detected in water samples from 69% of the tested shallow aquifer wells (Hayes et al., 1983). Another investigation by the U.S. Geological Survey in the southwestern tip of this region (McConnell et al., 1984) detected ethylene dibromide in the shallow aquifer. High levels of total nitrogen were also detected in groundwater quality samples from this area. Available information concerning the groundwater quality in the shallow aquifer clearly indicates that its water quality has been directly affected by the recent growth in agricultural activities. This pollution can cause a serious threat to the availability of good quality water which is supplied by an increasing reliance on groundwater. Between 1977 and 1980 the use of groundwater for irrigation increased by 62%. The main source of this water is the principal artesian aquifer, whose major recharge route includes the shallow aquifer. Establishment of a regional network for groundwater quality monitoring in the shallow aquifer is an important step in the protection and preservation of the principal artesian aquifer. This monitoring network can act as an early warning system for groundwater pollution problems in lower layers.
Hydrogeology of the shallow aquifer This section is primarily based on geohydrological information provided by Hayes et al. (1983). This report is the result of an extensive investigation of the hydrology of the principal artesian aquifer in the Dougherty Plain. This plain is bounded on the west by the Chattahoochee River, on the east by the Pelham Escarpment, and lies roughly southward of the updip limit of the principal artesian aquifer. It has a mild slope with an average land-surface elevation of about 160 ft above sea level. The Dougherty Plain is underlain by a succession of sand, clay, and carbonate rocks to a depth of more than 5000 ft. Its uppermost geologic units are: (1) the shallow aquifer (also known as the residuum); and (2) the principal artesian aquifer, as shown in Fig. 1. The shallow aquifer consists of a layer of sand and clay derived from chemical weathering of the Ocala Limestone. Clay content ranges from approximately 10 to 70%, with samples from 45 of 50 test wells consisting of more than 25% clay. The thickness of the shallow aquifer ranges from a few feet to slightly more than 125 ft, with an average thickness of approximately 50 ft. The shallow aquifer is not generally used as a source of water. It appears that its prime geohydrologic role is the recharge route to the principal artesian
91 aquifer. It is in this context that we study the desired network to monitor the quality of recharged water for the protection of the principal artesian aquifer. The shallow aquifer receives about 1.022 x 1012 gal of recharge per year, which occurs chiefly during J a n u a r y through May. Rainfall that is not evaporated, transpired, retained in the unsaturated zone as soil moisture, or discharged to streams, moves downward through the shallow aquifer to recharge the principal artesian aquifer. Preliminary results show that average daily recharge varies from about 0.1 to 2 Mgald lmi ~ Test drilling indicates that impermeable clay layers occur more commonly in the lower half of the shallow aquifer than in the upper half. On this basis, Hayes et al. (1983) suggest to consider the lower half of the shallow aquifer as a leaky layer. The thickness of this layer was combined with its vertical hydraulic conductivity to form the parameter known as leakance, as follows:
L(X) = Kv(X)/b(X)
(9)
where: L(X) = the leakance at X in ft d 1 ft 1 or l/d); Kv(X) the vertical hydraulic conductivity of the leaky (residuum confining) layer at X in ft d 1; and b(X) = the thickness of the leaky (residuum confining) layer at X in ft. Leakance is an indicator of potential recharge at any location from the shallow aquifer into the principal artesian aquifer. Estimated values of leakance vary from 0.00001 1/d to 0.36 l/d, with a median of 0.00172 l/d, as shown in Table 1. Leakance has been used as the main variable in the design of our monitoring network to identify potentially high recharge areas. STRUCTURAL ANALYSISAND MAPPINGRESULTS A structural analysis is conducted to determine the spatial statistics of leakance, which requires the selection of an appropriate distribution function for this variable. Freeze (1975) states that most field studies have indicated that the log-normal distribution is a suitable function to describe the statistical variations of the transmissivity data. The same argument has been used for hydraulic conductivity, too. Considering that the transmissivity is the product of the hydraulic conductivity and the saturated thickness, we can infer that the average saturated thickness is also log-normally distributed. This is based on the principle that the products or ratios of log-normally distributed variables are also log-normally distributed (Benjamin and Cornell, 1970). By extension, the leakance, which is the ratio of the vertical hydraulic conductivity and the thickness of the leaky layer (defined as the lower half of the shallow aquifer), is also log-normally distributed. This assumption is confirmed by our mapping results (Rouhani and Hall, 1987). Based on the above assumption, we can estimate different statistical properties of L(X) as: E[L(X)]
-
Var[L(X)]
exp{E[Y(X)] + Var[Y(X)]/2}
{E[L(X)]}~exp({Var[Y(X)]I- 1)
(10) (11)
92 TABLE 1 Estimated leakance for the shallow aquifer test wells (source Hayes et al., 1985) County
Well number
Well name
North (mi)
East* (mi)
Leakance (ft d ~ ft ~)
Baker
38 39 24 44 45 46 47 69 70 71 72 45 46 40 41 42 43 16 33 34 35 36 39 27 28 22 14 5 9
T. Renze RW JO-SU-LI TW1 B. Jordan TWl DPG J. Hall TW2 G. Bolton TW2 A. Newton School Bus Road TW1 Game and Fish TWl Nilo TW3 USMC Supply TW1 I. Newberry TW2 V. Evans Pied. Plant Farm TW1 S. Stocks TWl B. King TW1 H. Usry TW1 DP3 J. Fleet TW2 H. Meinders TW2 C. Bolton TW2 H. Davis TWl DP12 Roddenberry TW2 D. Harvey TW2 E. Stephens TW1 A. Vann TW1 DP9 C. Odom TW1
50.0 42.0 59.4 30.3 27.0 23.0 16.7 61.1 66.3 58.5 64.3 51.9 38.0 74.3 69.2 83.0 86.0 37.0 53.0 31.2 43.5 34.8 45.8 18.4 31.2 91.5 71.0 75.5 81.0
38.3 38.1 41.8 38.3 46.2 33.4 25.1 59.1 49.1 54.3 66.4 30.6 12.4 62.3 63.3 64.3 59.1 19.9 29.2 43.4 50.4 47.6 49.1 16.0 11.3 70.0 63.3 66.4 74.8
0.00762 0.00002 0.00003 0.00019 0.00010 0.00038 0.00755 0.00006 0.00042 0.00800 0.00007 0.00005 0.00013 0.00013 0.36000 0.00008 0.00001 0.00002 0.00049 0.00034 0.00010 0.00010 0.00003 0.00003 0.00007 0.00005 0.00001 0.00250 0.00003
Calhoun Decatur
Dougherty
Early Lee
Miller Mitchell
Seminole Sumpter Terrell Worth
* The origin corresponds to 30° 38'N, 85° 10'W.
m[L(X)]
=
exp {E[Y(X)]~
(12)
w h e r e : Y(X) - l o g - t r a n s f o r m o f L a t X ; a n d m = t h e m e d i a n . T h e r i s k v a l u e s o f L c a n a l s o b e c a l c u l a t e d b y s u b s t i t u t i n g L ( X ) f o r Z ( X ) i n e q n . (8). E s t i m a t i o n o f t h e o r d e r o f d r i f t (k) a s d e f i n e d i n e q n . (3), a n d , t h e c o v a r i a n c e parameters (C, al, a:~, a n d as), a s g i v e n i n e q n . (4), a r e c o n d u c t e d u s i n g a c o m p u t e r p r o g r a m b a s e d o n t h e w o r k s o f K a f r i t s a s a n d B r a s (1981). I t w a s determined that log-transformed leakance has a constant drift with a linear c o v a r i a n c e (k 0, K - 1.1552h). As a measure for the goodness-of-fit of the estimated covariance, we have u t i l i z e d a j a c k k n i f e e s t i m a t o r f o r ~: d e f i n e d a s t h e r a t i o o f t h e a c t u a l s u m o f s q u a r e d e r r o r s o f e s t i m a t i o n a n d t h e s u m o f k r i g i n g v a r i a n c e s ( R o u h a n i , 1983).
93
_
i
! L i
i
R
J
i
~rj LC o.
T
o.
o. o
o
f
R
0.0
]0. o
20.0
30,0
40.0
ERST
5C.0
50,0
7C.0
BO.O
90.0
(M ]LES }
Fig. 2. Map of expected leakance in ftd :ft
'
A perfect fit results in a v a l u e of 1 for iS. T h e c a l c u l a t e d v a l u e of 1.4822 displays a r e a s o n a b l e degree of goodness-of-fit, w h i c h is quite s a t i s f a c t o r y for field data.
Mapping results F o r m a p p i n g p u r p o s e s the D o u g h e r t y P l a i n a r e a is divided into a 20 x 22 grid with 5 mile i n c r e m e n t s . T h e a c t u a l m a p s are p r o d u c e d by the C o n t o u r p r o g r a m of D I S S P L A v e r s i o n 9.2 p a c k a g e , a v a i l a b l e at G e o r g i a T e c h ' s OCS Cyber computer. F i g u r e 2 displays the s p a t i a l d i s t r i b u t i o n s of the expected l e a k a n c e , w h i c h h a s a r e l a t i v e l y u n i f o r m v a l u e t h r o u g h o u t a l a r g e p o r t i o n of the s t u d y area. A sudden rise in the v a l u e of l e a k a n c e is i n d i c a t e d in the s o u t h e r n tip of this region. At this zone l e a k a n c e is a b o u t 4 to 10 times l a r g e r t h a n l e a k a n c e in o t h e r p a r t s of the plain. It m u s t be n o t e d t h a t due to the a s y m m e t r y of the l o g - n o r m a l d i s t r i b u t i o n some of the c a l c u l a t e d e x p e c t e d v a l u e s are excessively high. M e y e r (1975) s t a t e s that: " T h e e x p e c t a t i o n v a l u e is not so useful for an a s y m m e t r i c distribution. Often, m o r e significant are such m e a s u r e s as the median, mode, and the
94 C"
T F-Or"
3.0
10.3
2D.O
30.0
40.0
5C.0
50.0
70.C
80.0
90.0
EPST ~HiLES]
Fig. 3 Map of median leakance in ft d i ft geometric mean." So, we have utilized the median values, based on eqn. (12), as the representative map of this hydraulic variable, as shown in Fig. 3. In contrast to the expected value map (Fig. 2), this map displays moderate values with more spatial variations. The median map illustrates that the southern tip and the upper central zones are areas of high leakance. Figure 4 displays the spatial variations of estimation variances of the logtransformed leakance of the shallow aquifer, Var{ln[L] }. The boundary regions have higher variances due to the fact that these points are generally extrapolated, and thus, contain more uncertainty. As expected the southern and n o r t h er n tips show the highest levels of uncertainty. On the other hand, the middle of the upper half portion of the plain displays a lower variance of estimation, which is due to a relatively higher concentration of measurement points in this zone, see Fig. 6. Figure 5 illustrates the spatial distribution of the 10% risk values for the leakance (e.g., the chances that the actual values of leakance will be higher than these values are 10%). Comparison of Fig. 5 to the median map (Fig. 3) shows a clear increase in its values, with a lot less local variations. This smoothing is the result of combining estimated values with their variances.
95
o
i
/ i
z
o
J
o
0.0
lO.O
20.0
30.0
40.0
50.0
60.0
70.0
8J.O
90.0
ERST {MILES)
Fig. 4. Map of variance of log-transformed leakanee. SAMPLING SCHEMES Our proposed sampling schemes are based on three ranking criteria. The first one is a variance reduction criterion, which ranks the potential sampling points according to their variance reduction capabilities. This criterion is primarily concerned with the accuracy of the estimated values. The second criterion is the median ranking, which considers only the magnitude of the estimates. The third sampling criterion is the risk ranking, which is based on both the accuracy and the magnitude of the estimated values. For sampling purposes we have defined the 32 points shown in Fig. 6 as the potential sampling sites. These points are distributed uniformly over the Dougherty Plain. This figure also shows the location of existing sampling sites which are scattered t h r o u g h o u t this region. There are a few concentrations of data points, such as in the middle of the upper portion. However in general, their distribution can be considered relatively uniform. Such a distribution allows us to examine the relative importance of potential sampling sites without any implicit bias. The sampling activities are conducted in a nonsequential manner, which implies the constancy of the covariance function.
96
Ld
~J Z
O.@
]0.9
20.0
30.0
40.0
EFST
50.U
BO.U
/U.U
oJ.u
uu.C
(NILES~
Fig. 5. Map of 10°4 risk of leakance in ft d 1ft l This assumption can be violated if the m e a s u r e d values at the new sampling sites t u r n out to be significantly different from t h e i r predicted values. F o r a study of the resilience of v a r i a n c e r e d u c t i o n analysis u n d e r sequential samplings readers are referred to R o u h a n i and Fiering (1986). At each r o u n d of kriging the site a m o n g p o t e n t i a l locations with the maximal v a l u e of the stated criterion is selected as the next m e a s u r e m e n t location. This yields a sequence of n points a m o n g the 32 p o t e n t i a l sites for f u r t h e r sampling. As noted in R o u h a n i and Fiering (1986), the r a n k list can be used as a shopping list; we utilize it, from the top down, until the budget is e x h a u s t e d or some o t h e r condition is met.
Sampling based on the variance reduction analysis F i g u r e 7 shows the sequence of selected points based on v a r i a n c e r e d u c t i o n analysis. This sequence is a valid sampling guideline as long as our sole objective is the m a x i m i z a t i o n of i n c r e m e n t a l information. As seen, the nodes at the e a s t e r n and s o u t h e r n b o u n d a r i e s have h i g h e r ranks, which implies t h a t
97
100
.
90
O°O •
80
0 0 • (~•0 ° (~0•© 0 ( ~ 0 0
70 6O E "r 5 0 t0 z
~O
4O
0 0 ~)'• 0 o .o.o • o O0
30 20 10 0
o
1'o
2'o 3'0 4'0
5'0
6'0 7'0 8'o 90
E A S T (miles)
Fig. 6. Potential and existing sampling points (solid circles are existing data points). 100
@@ ® @ @®@® @®@® @@®@ ®@@@@ @®®®® @@®® @Q
90 80 70
60"~ "l- 50 nO 4O z 30 20 10
0
1'0
2'0 3'0 4'o
5'0
6'0 7'o
EAST (miles) Fig. 7. Variance reduction analysis for leakance.
go
9o
98
200-
(l/dr
OTV
100
)
1'0
2'0
3'0
RANK Fig. 8. TOTV and TVR at each round of kriging for leakance. higher priorities must be assigned to these zones when planning a monitoring network. The central zone and the western boundary, on the other hand, have lower ranks which means that they should not be considered as primary sampling locations. Figure 8 illustrates the level of u n c e r t a i n t y at each round of sampling in terms of T O T V (i.e., total sum of estimation variances), and the amount of information gain due to each new sampling in terms of its T V R , eqn. (7). As shown, sampling at the top five points yields the highest amount of information gain, while additional measurements appear to produce only marginal benefits. One can assume that there must be a finite number of new sampling sites, such that any additional sampling will result in information gains that cannot be justified in economic terms. S a m p l i n g based on median r a n k i n g
Median ranking, shown in Fig. 9, provides an alternative criterion based on the estimated magnitudes of leakance. However, the accuracy of the estimated
99 100
70
@@ ® @ ®@@@
60
@@@@
90 80
E "~ l- 50 nO 40 z
@@®® ®@@®@ @@@®® ®®@Q
30 20
®@
10 0
1'0
20
3'0
4'0
5'0
6'0
7'0
8'0
90
EAST (miles) Fig. 9. Median ranking for leakance.
100"
90
® ®
80
®
70
@®®@
~6o
@@@®
l.-rr
®@@@@
040 z
30
® ® @ @ ®
20
®@@Q
,o
•
@
20
3'0
0
@
0
1'0
4'0
5'0
EAST (miles) Fig. 10. 10% risk ranking for leakance.
6'0
7'0
8'0
90
1O0
100
O@ O0 @0®@ 0@00 000(~ 00000 Q 0000 00@@ QQ
90 80 70 60 E v 50 -r knO 40 Z 30 20 10 0
1'0
2'0
3'0
4'0
5'0
6'0
7~0
8'0
90
EAST (miles) Fig. 11. 1%, risk ranking for leakance.
100-
O0 0 0 O@~)Q 0@00 0@00 00000 ®000© 00@@ @@
90 80 70 60 ~
50
0 z
40
I
30 20 10
0
10
2'o 3'0 4'0
s'o 0'0 7'0 8'0 90
EAST (miles) Fig. 12.40%, risk ranking for leakance.
101
field has been ignored. In this ranking the southern tip and the upper central portions of the Dougherty Plain have the highest priorities, while the lower central and the n o r t h e r n tip zones and the eastern boundary have acquired lower priorities. Sampling based on risk ranking Risk value which depends on both the magnitude and the accuracy of an estimated leakance, seems to be a more suitable criterion for our sampling purposes. As illustrated in Fig. 10, the southern tip and the central zones have been identified as the primary candidates for sampling. The risk ranking, unlike the variance reduction analysis, has not ignored areas of high leakage with good accuracy, such as the central region (nodes ranked 4th, 5th, and 6th). On the other hand, areas with low potential recharge, such as the lower central and the n o r th er n tip zones along with the western boundary, have received lower priorities. These results demonstrate the advantages of risk ranking, which in addition to the variance, considers the actual magnitudes of leakance, and yields a more balance guideline for sampling. Equation (8) shows that the risk value is basically a weighted sum of the expected value and the estimation variance of the log-transformed leakance. Thus, as we decrease the probability of exceedence of our risk values, we are in fact giving more weight to the variance, and conversely. For instance, the 1% risk ranking, shown in Fig. 11, displays a move toward the results of variance reduction analysis as shown in Fig. 7. On the other hand, as expected, the result of the 40% risk ranking, shown in Fig. 12, is almost identical to the median ranking given in Fig. 9. The risk ranking gives us a degree of flexibility to adjust this weighted sum according to our specific sampling objectives. This work clearly demonstrates the potentials of geostatistical techniques in the design of regional groundwater monitoring networks. Following these optimally oriented procedures can lead to more efficient sampling schemes, at lower costs. ACKNOWLEDGMENTS This research was supported by grant ECE-8503897 from the National Science Foundation, and U.S. Department of Interior, Geological Survey Project G-1219(05). The authors would like to thank Prof. Ram Arora of Georgia State University for his help and guidance.
REFERENCES Benjamin, J.R. and Cornell, C.A., 1970. Probability, Statistics, and Decision for Civil Engineers. McGraw-Hill, New York, N.Y. Delfiner, P., 1975. Linear estimation of non-stationary spatial phenomena. In: M. Guarascio. M. David, and C. Huijbregts (Editors), AdvancedGeostatistics in the Mining Industry. Nato Adv. Stud. Inst.. Reidel, Hingham, Mass.
102 Freeze, R.A., 1975. A stochastic-conceptual analysis of one-dimensional groundwater flow in nonuniform homogeneous media. Water Resour. Res., 11(5): 72~741. Georgia Department of Agriculture, 1985. Georgia Agricultural Facts. Ga. Dep. Agric., Atlanta, Ga. Hayes, L.R., Maslia, M.L. and Meeks, W.C., 1983. Hydrology and model evaluation of the principal artesian aquifer, Dougherty Plain, Southwest Georgia. Bull. 97, Dep. Nat. Resour.. Environ. Prot. Div., Ga. Geol. Surv., Atlanta, Ga. Journel, A.G. and Huijbregts, C., 1978. Mining Geostatistics. Academic Press, London, 600 pp. Kaf?itsas, J. and Bras, R.L., 1981. The practice of kriging. Rep. No. 263, Ralph M. Parsons I,ab.. M.I.T., Cambridge, Mass. McConnel, J.B., Hicks, D.W., Lowe, L.E., Cohen, S.Z. and Jovanovich, A.P., 1984. Investigation of ethylene dibromide (EDB) in ground water in Seminole County, Georgia. U.S. Dep. Inter., Geol. Surv., Circ. 933, Reston, Va. Meyer, S.T., 1975. Data Analysis for Scientists and Engineers. Wiley, New York, N.Y., 285 pp. Rouhani, S.. 1983. Optimal data collection in random fields. Ph.D. Thesis, Div. Appl. Sci., Harvard Univ., Cambridge, Mass. Rouhani, S., 1985. Variance reduction analysis. Water Resour. Res., 21(6): 837-846. Rouhani, S. and Fiering, M.B., 1986. Resilience of a statistical sampling scheme. J. Hydrol., 89: 1 11. Rouhani, S. and Hall, T.J., 1987. Optimal schemes for ground water quality monitoring in the shallow aquifer, Dougherty Plain, southwestern Georgia. Tech. Complet. Rep., USDI/USGS Proj. G-1219 (05), ERC 03 87, Environ. Resour. Center, Georgia Inst. Technol., Atlanta, Ga., 49 pp.
85 Printed in The Netherlands
[4] G E O S T A T I S T I C A L S C H E M E S FOR G R O U N D W A T E R S A M P L I N G
SHAHROKH ROUHANI and TIMOTHY J. HALL
School of Civil Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332 (U.S.A.) (Received March 12, 1987; revised and accepted January 8, 1988)
ABSTRACT Rouhani, S. and Hall, T.J., 1988. Geostatistical schemes for groundwater sampling. J. Hydrol., 103:85 102. Geostatistical techniques offer efficient tools for design of groundwater sampling networks. They include procedures for the selection of the best sequences of sampling points, such as: variance reduction analysis, median ranking, and risk ranking. Variance reduction analysis considers primarily the accuracy of the estimated field, while median ranking is based only on the magnitude of the estimated values. Risk ranking is a compromise between these procedures that appears to yield a more balanced guideline for cases when planners desire to acquire maximum information, while monitoring areas where the variable of interest exhibits critical values. These procedures are used for the design of a regional shallow groundwater quality monitoring network in the Dougherty Plain, located in southwest Georgia. The shallow aquifer of concern is the main recharge route to a semiconfined aquifer which is the primary source of water in this region. The desired monitoring network acts as an early warning system for groundwater po]lution in deeper layers. Leakance data are utilized to identify the primary sampling locations. Statistical analyses indicate that leakance has a log-normal distribution with a constant drift and a linear spatial covariance. The results of our risk rankings demonstrate that the southern tip of the Dougherty Plain and its upper central zone should be the prime targets of our monitoring activities.
INTRODUCTION
An increasing number of researchers are advocating the use of geostatistical techniques in the design of groundwater sampling networks (Rouhani and Fiering, 1986). For this purpose the sampled variables should be viewed as stochastic processes. In particular, spatial variables such as the transmissivity, storativity, and steady-state piezometric heads may be considered as random fields, with the following structure: Z(X,)
=
M(X~) + R ( ~ )
(1)
where: Z(X/) = the actual value of the variable of interest at X/; M(X/) = a slowly varying deterministic function known as the drift which is equal to the expected value of Z at X~; and R(X~) = a spatially fluctuating random function with zero expectation and a spatial covariance function. As Rouhani (1985) states, the data management of geohydrologic random
0022-1694/88/$03.50
© 1988 Elsevier Science Publishers B.V.
86 fields depends upon the objectives of the planning approach. Often very little data are available. Furthermore, the measured values may be clustered together and therefore not provide information about the whole field. In such situations, planners may wish to design a data collection scheme in order to better define their variables of interest. The questions that immediately arise are: How should we design such a regional sampling network? More specifically, what criteria should be utilized as the basis of our network design? Where are the best locations for sampling sites? Do the selected sites remain unchanged if we change our sampling criterion? Geostatistical techniques such as kriging provide tools to answer these questions. Commonly geostatistical sampling procedures are based on the maximization of incremental information subject to budget constraints. For instance, Rouhani (1985) proposes variance reduction analysis to select those sequences of n points so chosen from m potential sites to maximize reduction in the total variance of estimates. In these procedures the criteria of site selection are such that they give more priority to points with high estimation variances, regardless of their estimated magnitudes. These criteria are suitable for cases where the estimated value of the variable of concern is not of primary importance. In many cases, however, planners are not only interested in gaining as much information as possible, but also to monitor areas where the sampled variable exhibits critical values. This requires some modifications in the selection criterion to include both the accuracy and the magnitude of the estimated values. To accomplish this task we propose to utilize the risk value as an alternative selection criterion, as described in the following sections. GEOSTATIST1CALANALYSIS Geostatistical techniques provide means for the statistical analysis of spatially distributed random variables. Matheron (Journel and Huijbregts, 1978) states that: ~Geostatistics is the application of the formalism of random functions to the reconnaissance and estimation of natural phenomena." At first glance, the measured values of many geohydrological variables appear to be random and unpredictable. However, values taken at neighboring points may be related by statistical functions. Geostatistics enables us to estimate these statistical relations and utilize them to estimate the value of the variable at an arbitrary point. Kriging is one of the important geostatistical techniques which has been applied to many geohydrologic problems for estimation and sampling purposes. This algorithm is basically a statistical interpolation technique. In point kriging one estimates the value of a random field at any arbitrary point. X~), from values measured at other points, in a linear form of: N
Z(Xo) = ~ )~ioZ(X~) i
1
(2)
87
where: Z(X0) = k r i g i n g e s t i m a t e at 3/o; Z(Xi) - m e a s u r e d v a l u e at Xi, i - 1 ..... N; 2li0 - k r i g i n g w e i g h t for Z(X~) to e s t i m a t e Z(X0); a n d N - n u m b e r of measu r e m e n t points. T h e 2i0 are defined by two criteria: (1) unbiasedness; E[2(X0) Z(X~,)] - 0, w h e r e Z(Xo) is the a c t u a l v a l u e of the field at X0; and (2) m i n i m u m s q u a r e d e r r o r of estimation: this r e q u i r e s E[Z(X0) - Z(X0)] 2 to be m i n i m u m . The m i n i m u m s q u a r e d e r r o r of e s t i m a t i o n is also a m e a s u r e for the a c c u r a c y of o u r e s t i m a t e s which is k n o w n as e s t i m a t i o n (or kriging) v a r i a n c e . K r i g i n g a s s u m e s t h a t the v a r i a b l e of i n t e r e s t is a r a n d o m field w i t h a s t r u c t u r e s i m i l a r to eqn. (1). In this w o r k w i t h o u t a n y loss of g e n e r a l i t y , it is a s s u m e d t h a t the e x p e c t e d v a l u e of Z(X0), the so-called drift, is a p o l y n o m i a l of kth order: M(Xo)
l(k) ~ bJp(Xo)
, p
(3)
1
where: b t , - fixed u n k n o w n coefficients; f p ( X ) = the p t h m o n o m i a l ; and l(k) = the n u m b e r of m o n o m i a l s in the a b o v e kth o r d e r polynomial. In a t w o - d i m e n s i o n a l s p a c e with C a r t e s i a n c o o r d i n a t e s (x,y), a s e c o n d - o r d e r p o l y n o m i a l (k - 2, and 1 - 6) h a s the following m o n o m i a l s : 1, x, y, xy, x 2, and y2.
T h e r a n d o m c o m p o n e n t in eqn. (1), R, is a s s u m e d to h a v e a zero m e a n and a c o v a r i a n c e , defined as K[R(Xi),R(X~)] = E [ R ( X i ) R ( X j ) ]. It is f u r t h e r a s s u m e d t h a t the c o v a r i a n c e is an isotropic f u n c t i o n of the d i s t a n c e b e t w e e n a n y two points: K 0 = K([X/ - Xj[). T h e e s t i m a t i o n of the c o v a r i a n c e is a c c o m p l i s h e d u s i n g the s t r u c t u r a l a n a l y s i s p r o p o s e d by Delfiner (1975). The s u g g e s t e d p o l y n o m i a l c o v a r i a n c e f u n c t i o n h a s the following form: k
K(h)
Cb(h) +
~ p
av.
,h 2p'1
(4)
1
where: K - the c o v a r i a n c e function; h = the length of v e c t o r d i s t a n c e b e t w e e n a n y two points; ap - the p t h coefficient; C = the point v a r i a n c e ; = the D i r a c d e l t a function; and k = the o r d e r of the p o l y n o m i a l drift. T h e e s t i m a t i o n v a r i a n c e c a n t h e n be c a l c u l a t e d as:
Po i
N
N
2
2 0j
:
)'io~tjoK(IXi - XJt )
(5)
0
w h e r e Vo is the e s t i m a t i o n v a r i a n c e at X0, Zoo - - 1, and the r e s t of definitions are a l r e a d y given. F o r f u r t h e r i n f o r m a t i o n is r e f e r r e d to R o u h a n i (1983 and 1985). V a r i a n c e reduction a n a l y s i s T h e e s t i m a t i o n v a r i a n c e in eqn. (5) is a m e a s u r e for the a c c u r a c y of the e s t i m a t e d v a l u e at X0. Utilizing this function, v a r i a n c e r e d u c t i o n a n a l y s i s provides a s c h e m e based on the m a x i m i z a t i o n of i n c r e m e n t a l i n f o r m a t i o n
88 subject to b u d g e t c o n s t r a i n t s (Rouhani, 1985). This p r o c e d u r e first m e a s u r e s the m a g n i t u d e of i n f o r m a t i o n g a i n at a n y e s t i m a t e d point, X0, due to a m e a s u r e m e n t at X * (the a r b i t r a r y l o c a t i o n of a p o t e n t i a l s a m p l i n g site), in t e r m s of r e d u c t i o n s in V0 a f t e r to the new sampling, as: N
VR,~.
ltk)
[Ko. = ~ i
2i, K,o1
~ p
l,p.fp(Xo)['/[V.(N)]
(6)
1
where: VR.+ = the v a r i a n c e r e d u c t i o n (i.e., a c c u r a c y or i n f o r m a t i o n gain) at X¢~ due to a m e a s u r e m e n t at X,; Ki~ the c o v a r i a n c e f u n c t i o n b e t w e e n X, and Xj; 2,+ = the o p t i m a l w e i g h t of Z(Xi) in e s t i m a t i o n of Z(X.) p r i o r to the new m e a s u r e m e n t ; /tp, = the L a g r a n g e m u l t i p l i e r for the p t h m o n o m i a l filtering c o n s t r a i n t in the k r i g i n g s y s t e m for the e s t i m a t i o n of Z(X.) p r i o r to the new m e a s u r e m e n t ; N = the n u m b e r of existing d a t a points prior to the new meas u r e m e n t ; and V,(N) the k r i g i n g v a r i a n c e at X, p r i o r to the new m e a s u r e ment. It t h e n e x p a n d s the a b o v e m e a s u r e to c o v e r the e n t i r e field as follows:
TVR.
~ VRj,
(7)
i
where: TVR, = the t o t a l v a r i a n c e r e d u c t i o n due to a m e a s u r e m e n t at X,; and j - the set of e s t i m a t e d points. TVR, r e p r e s e n t s the a m o u n t of r e g i o n a l inf o r m a t i o n g a i n due to a new m e a s u r e m e n t at X,. At e a c h r o u n d of sampling, the site a m o n g p o t e n t i a l s a m p l i n g l o c a t i o n s with m a x i m a l TVR, is selected as the n e x t m e a s u r e m e n t location. This yields a s e q u e n c e of n points a m o n g m sites for f u r t h e r sampling.
Risk ranking T h e a b o v e a p p r o a c h gives m o r e p r i o r i t y to points with h i g h e s t i m a t i o n v a r i a n c e s , r e g a r d l e s s of t h e i r e s t i m a t e d m a g n i t u d e . S u c h a c r i t e r i o n is not s u i t a b l e for a case w h e n p l a n n e r s desire to g a i n m a x i m u m i n f o r m a t i o n , while m o n i t o r i n g a r e a s w h e r e the v a r i a b l e of c o n c e r n exhibits critical values. F o r i n s t a n c e , in o u r case study, the desired n e t w o r k is aimed to m o n i t o r the q u a l i t y of r e c h a r g e d w a t e r at critical l e a k y points. So, we are not o n l y i n t e r e s t e d in g a i n i n g as m u c h i n f o r m a t i o n as possible, b u t also to m o n i t o r a r e a s with p o t e n t i a l l y h i g h levels of r e c h a r g e . This implies t h a t the s e l e c t i o n c r i t e r i o n should be modified to include b o t h the a c c u r a c y a n d the m a g n i t u d e of o u r e s t i m a t e d values. T h i s o b j e c t i v e c a n be a c c o m p l i s h e d by utilizing the risk v a l u e as an altern a t i v e selection criterion. T h e risk v a l u e is defined as the v a l u e of the v a r i a b l e of i n t e r e s t w h o s e p r o b a b i l i t y of e x c e e d e n c e is ~ percent. F o r a l o g - n o r m a l l y d i s t r i b u t e d v a r i a b l e the risk v a l u e can be w r i t t e n as: Z~(X) =
exp[E{ln[Z(X)]} + z, (Var{ln[Z(X)l}) 1'~]
(8)
where: Z~(X) = the risk v a l u e at X w i t h the p r o b a b i l i t y of e x c e e d e n c e of a
89
.
ANTA
DOUGHERTY
Fig. 1. Geological profile of Georgia, U.S.A. percent; Z(X) = r a n d o m variable of interest at X; Var( ) = the variance; and z~ - the standardized n o r m a l l y distributed r a n d o m variable with a probability of e x c e e d e n c e of ~-percent. In o u r case study, at each r o u n d of kriging, the point with maximal risk value is identified as the n e x t best sampling site, which will yield an ordered list of points for site selections. As we r e d u c e the level of o u r risk, ~, more weight will be given to the variance. It can be inferred t h a t v a r i a n c e r e d u c t i o n analysis is an e x t r e m e case of risk r a n k i n g when v a r i a n c e has the complete weight. On the o t h e r extreme, when ~ = 50%, we h a v e the m e d i a n r a n k i n g which totally ignores the a c c u r a c y of the estimates. APPLICATION For o u r case study we have selected the g r o u n d w a t e r q u a l i t y m o n i t o r i n g in the shallow aquifer of the D o u g h e r t y Plain in s o u t h w e s t e r n Georgia, as shown in Fig. 1. This a r e a is a rapidly growing a g r i c u l t u r a l region, which is u n d e r l a i n by a succession of sand, clay and c a r b o n a t e rocks to a depth of more t h a n 5000 ft, forming one of the most p r o d u c t i v e m u l t i l a y e r aquifer systems in the United States. The a b u n d a n t supply of good quality g r o u n d w a t e r , along with a mild climate, and a flat to gently rolling t e r r a i n has induced significant g r o w t h in the a g r i c u l t u r a l p r o d u c t i o n of this region. This growth is accom-
90 panied by an increased use of fertilizers and pesticides. Recent data indicates that the average fertilizer use in the study area is approximately 105 t acre In 1977 this level was even higher at about 150 t acre l (Georgia Department of Agriculture, 1985). Some of the chemical components of these fertilizers and pesticides are toxic to human, long-lasting, and tend to accumulate in the hydrogeological system. In a recent study by the Georgia Department of Natural Resources, pesticides were detected in water samples from 69% of the tested shallow aquifer wells (Hayes et al., 1983). Another investigation by the U.S. Geological Survey in the southwestern tip of this region (McConnell et al., 1984) detected ethylene dibromide in the shallow aquifer. High levels of total nitrogen were also detected in groundwater quality samples from this area. Available information concerning the groundwater quality in the shallow aquifer clearly indicates that its water quality has been directly affected by the recent growth in agricultural activities. This pollution can cause a serious threat to the availability of good quality water which is supplied by an increasing reliance on groundwater. Between 1977 and 1980 the use of groundwater for irrigation increased by 62%. The main source of this water is the principal artesian aquifer, whose major recharge route includes the shallow aquifer. Establishment of a regional network for groundwater quality monitoring in the shallow aquifer is an important step in the protection and preservation of the principal artesian aquifer. This monitoring network can act as an early warning system for groundwater pollution problems in lower layers.
Hydrogeology of the shallow aquifer This section is primarily based on geohydrological information provided by Hayes et al. (1983). This report is the result of an extensive investigation of the hydrology of the principal artesian aquifer in the Dougherty Plain. This plain is bounded on the west by the Chattahoochee River, on the east by the Pelham Escarpment, and lies roughly southward of the updip limit of the principal artesian aquifer. It has a mild slope with an average land-surface elevation of about 160 ft above sea level. The Dougherty Plain is underlain by a succession of sand, clay, and carbonate rocks to a depth of more than 5000 ft. Its uppermost geologic units are: (1) the shallow aquifer (also known as the residuum); and (2) the principal artesian aquifer, as shown in Fig. 1. The shallow aquifer consists of a layer of sand and clay derived from chemical weathering of the Ocala Limestone. Clay content ranges from approximately 10 to 70%, with samples from 45 of 50 test wells consisting of more than 25% clay. The thickness of the shallow aquifer ranges from a few feet to slightly more than 125 ft, with an average thickness of approximately 50 ft. The shallow aquifer is not generally used as a source of water. It appears that its prime geohydrologic role is the recharge route to the principal artesian
91 aquifer. It is in this context that we study the desired network to monitor the quality of recharged water for the protection of the principal artesian aquifer. The shallow aquifer receives about 1.022 x 1012 gal of recharge per year, which occurs chiefly during J a n u a r y through May. Rainfall that is not evaporated, transpired, retained in the unsaturated zone as soil moisture, or discharged to streams, moves downward through the shallow aquifer to recharge the principal artesian aquifer. Preliminary results show that average daily recharge varies from about 0.1 to 2 Mgald lmi ~ Test drilling indicates that impermeable clay layers occur more commonly in the lower half of the shallow aquifer than in the upper half. On this basis, Hayes et al. (1983) suggest to consider the lower half of the shallow aquifer as a leaky layer. The thickness of this layer was combined with its vertical hydraulic conductivity to form the parameter known as leakance, as follows:
L(X) = Kv(X)/b(X)
(9)
where: L(X) = the leakance at X in ft d 1 ft 1 or l/d); Kv(X) the vertical hydraulic conductivity of the leaky (residuum confining) layer at X in ft d 1; and b(X) = the thickness of the leaky (residuum confining) layer at X in ft. Leakance is an indicator of potential recharge at any location from the shallow aquifer into the principal artesian aquifer. Estimated values of leakance vary from 0.00001 1/d to 0.36 l/d, with a median of 0.00172 l/d, as shown in Table 1. Leakance has been used as the main variable in the design of our monitoring network to identify potentially high recharge areas. STRUCTURAL ANALYSISAND MAPPINGRESULTS A structural analysis is conducted to determine the spatial statistics of leakance, which requires the selection of an appropriate distribution function for this variable. Freeze (1975) states that most field studies have indicated that the log-normal distribution is a suitable function to describe the statistical variations of the transmissivity data. The same argument has been used for hydraulic conductivity, too. Considering that the transmissivity is the product of the hydraulic conductivity and the saturated thickness, we can infer that the average saturated thickness is also log-normally distributed. This is based on the principle that the products or ratios of log-normally distributed variables are also log-normally distributed (Benjamin and Cornell, 1970). By extension, the leakance, which is the ratio of the vertical hydraulic conductivity and the thickness of the leaky layer (defined as the lower half of the shallow aquifer), is also log-normally distributed. This assumption is confirmed by our mapping results (Rouhani and Hall, 1987). Based on the above assumption, we can estimate different statistical properties of L(X) as: E[L(X)]
-
Var[L(X)]
exp{E[Y(X)] + Var[Y(X)]/2}
{E[L(X)]}~exp({Var[Y(X)]I- 1)
(10) (11)
92 TABLE 1 Estimated leakance for the shallow aquifer test wells (source Hayes et al., 1985) County
Well number
Well name
North (mi)
East* (mi)
Leakance (ft d ~ ft ~)
Baker
38 39 24 44 45 46 47 69 70 71 72 45 46 40 41 42 43 16 33 34 35 36 39 27 28 22 14 5 9
T. Renze RW JO-SU-LI TW1 B. Jordan TWl DPG J. Hall TW2 G. Bolton TW2 A. Newton School Bus Road TW1 Game and Fish TWl Nilo TW3 USMC Supply TW1 I. Newberry TW2 V. Evans Pied. Plant Farm TW1 S. Stocks TWl B. King TW1 H. Usry TW1 DP3 J. Fleet TW2 H. Meinders TW2 C. Bolton TW2 H. Davis TWl DP12 Roddenberry TW2 D. Harvey TW2 E. Stephens TW1 A. Vann TW1 DP9 C. Odom TW1
50.0 42.0 59.4 30.3 27.0 23.0 16.7 61.1 66.3 58.5 64.3 51.9 38.0 74.3 69.2 83.0 86.0 37.0 53.0 31.2 43.5 34.8 45.8 18.4 31.2 91.5 71.0 75.5 81.0
38.3 38.1 41.8 38.3 46.2 33.4 25.1 59.1 49.1 54.3 66.4 30.6 12.4 62.3 63.3 64.3 59.1 19.9 29.2 43.4 50.4 47.6 49.1 16.0 11.3 70.0 63.3 66.4 74.8
0.00762 0.00002 0.00003 0.00019 0.00010 0.00038 0.00755 0.00006 0.00042 0.00800 0.00007 0.00005 0.00013 0.00013 0.36000 0.00008 0.00001 0.00002 0.00049 0.00034 0.00010 0.00010 0.00003 0.00003 0.00007 0.00005 0.00001 0.00250 0.00003
Calhoun Decatur
Dougherty
Early Lee
Miller Mitchell
Seminole Sumpter Terrell Worth
* The origin corresponds to 30° 38'N, 85° 10'W.
m[L(X)]
=
exp {E[Y(X)]~
(12)
w h e r e : Y(X) - l o g - t r a n s f o r m o f L a t X ; a n d m = t h e m e d i a n . T h e r i s k v a l u e s o f L c a n a l s o b e c a l c u l a t e d b y s u b s t i t u t i n g L ( X ) f o r Z ( X ) i n e q n . (8). E s t i m a t i o n o f t h e o r d e r o f d r i f t (k) a s d e f i n e d i n e q n . (3), a n d , t h e c o v a r i a n c e parameters (C, al, a:~, a n d as), a s g i v e n i n e q n . (4), a r e c o n d u c t e d u s i n g a c o m p u t e r p r o g r a m b a s e d o n t h e w o r k s o f K a f r i t s a s a n d B r a s (1981). I t w a s determined that log-transformed leakance has a constant drift with a linear c o v a r i a n c e (k 0, K - 1.1552h). As a measure for the goodness-of-fit of the estimated covariance, we have u t i l i z e d a j a c k k n i f e e s t i m a t o r f o r ~: d e f i n e d a s t h e r a t i o o f t h e a c t u a l s u m o f s q u a r e d e r r o r s o f e s t i m a t i o n a n d t h e s u m o f k r i g i n g v a r i a n c e s ( R o u h a n i , 1983).
93
_
i
! L i
i
R
J
i
~rj LC o.
T
o.
o. o
o
f
R
0.0
]0. o
20.0
30,0
40.0
ERST
5C.0
50,0
7C.0
BO.O
90.0
(M ]LES }
Fig. 2. Map of expected leakance in ftd :ft
'
A perfect fit results in a v a l u e of 1 for iS. T h e c a l c u l a t e d v a l u e of 1.4822 displays a r e a s o n a b l e degree of goodness-of-fit, w h i c h is quite s a t i s f a c t o r y for field data.
Mapping results F o r m a p p i n g p u r p o s e s the D o u g h e r t y P l a i n a r e a is divided into a 20 x 22 grid with 5 mile i n c r e m e n t s . T h e a c t u a l m a p s are p r o d u c e d by the C o n t o u r p r o g r a m of D I S S P L A v e r s i o n 9.2 p a c k a g e , a v a i l a b l e at G e o r g i a T e c h ' s OCS Cyber computer. F i g u r e 2 displays the s p a t i a l d i s t r i b u t i o n s of the expected l e a k a n c e , w h i c h h a s a r e l a t i v e l y u n i f o r m v a l u e t h r o u g h o u t a l a r g e p o r t i o n of the s t u d y area. A sudden rise in the v a l u e of l e a k a n c e is i n d i c a t e d in the s o u t h e r n tip of this region. At this zone l e a k a n c e is a b o u t 4 to 10 times l a r g e r t h a n l e a k a n c e in o t h e r p a r t s of the plain. It m u s t be n o t e d t h a t due to the a s y m m e t r y of the l o g - n o r m a l d i s t r i b u t i o n some of the c a l c u l a t e d e x p e c t e d v a l u e s are excessively high. M e y e r (1975) s t a t e s that: " T h e e x p e c t a t i o n v a l u e is not so useful for an a s y m m e t r i c distribution. Often, m o r e significant are such m e a s u r e s as the median, mode, and the
94 C"
T F-Or"
3.0
10.3
2D.O
30.0
40.0
5C.0
50.0
70.C
80.0
90.0
EPST ~HiLES]
Fig. 3 Map of median leakance in ft d i ft geometric mean." So, we have utilized the median values, based on eqn. (12), as the representative map of this hydraulic variable, as shown in Fig. 3. In contrast to the expected value map (Fig. 2), this map displays moderate values with more spatial variations. The median map illustrates that the southern tip and the upper central zones are areas of high leakance. Figure 4 displays the spatial variations of estimation variances of the logtransformed leakance of the shallow aquifer, Var{ln[L] }. The boundary regions have higher variances due to the fact that these points are generally extrapolated, and thus, contain more uncertainty. As expected the southern and n o r t h er n tips show the highest levels of uncertainty. On the other hand, the middle of the upper half portion of the plain displays a lower variance of estimation, which is due to a relatively higher concentration of measurement points in this zone, see Fig. 6. Figure 5 illustrates the spatial distribution of the 10% risk values for the leakance (e.g., the chances that the actual values of leakance will be higher than these values are 10%). Comparison of Fig. 5 to the median map (Fig. 3) shows a clear increase in its values, with a lot less local variations. This smoothing is the result of combining estimated values with their variances.
95
o
i
/ i
z
o
J
o
0.0
lO.O
20.0
30.0
40.0
50.0
60.0
70.0
8J.O
90.0
ERST {MILES)
Fig. 4. Map of variance of log-transformed leakanee. SAMPLING SCHEMES Our proposed sampling schemes are based on three ranking criteria. The first one is a variance reduction criterion, which ranks the potential sampling points according to their variance reduction capabilities. This criterion is primarily concerned with the accuracy of the estimated values. The second criterion is the median ranking, which considers only the magnitude of the estimates. The third sampling criterion is the risk ranking, which is based on both the accuracy and the magnitude of the estimated values. For sampling purposes we have defined the 32 points shown in Fig. 6 as the potential sampling sites. These points are distributed uniformly over the Dougherty Plain. This figure also shows the location of existing sampling sites which are scattered t h r o u g h o u t this region. There are a few concentrations of data points, such as in the middle of the upper portion. However in general, their distribution can be considered relatively uniform. Such a distribution allows us to examine the relative importance of potential sampling sites without any implicit bias. The sampling activities are conducted in a nonsequential manner, which implies the constancy of the covariance function.
96
Ld
~J Z
O.@
]0.9
20.0
30.0
40.0
EFST
50.U
BO.U
/U.U
oJ.u
uu.C
(NILES~
Fig. 5. Map of 10°4 risk of leakance in ft d 1ft l This assumption can be violated if the m e a s u r e d values at the new sampling sites t u r n out to be significantly different from t h e i r predicted values. F o r a study of the resilience of v a r i a n c e r e d u c t i o n analysis u n d e r sequential samplings readers are referred to R o u h a n i and Fiering (1986). At each r o u n d of kriging the site a m o n g p o t e n t i a l locations with the maximal v a l u e of the stated criterion is selected as the next m e a s u r e m e n t location. This yields a sequence of n points a m o n g the 32 p o t e n t i a l sites for f u r t h e r sampling. As noted in R o u h a n i and Fiering (1986), the r a n k list can be used as a shopping list; we utilize it, from the top down, until the budget is e x h a u s t e d or some o t h e r condition is met.
Sampling based on the variance reduction analysis F i g u r e 7 shows the sequence of selected points based on v a r i a n c e r e d u c t i o n analysis. This sequence is a valid sampling guideline as long as our sole objective is the m a x i m i z a t i o n of i n c r e m e n t a l information. As seen, the nodes at the e a s t e r n and s o u t h e r n b o u n d a r i e s have h i g h e r ranks, which implies t h a t
97
100
.
90
O°O •
80
0 0 • (~•0 ° (~0•© 0 ( ~ 0 0
70 6O E "r 5 0 t0 z
~O
4O
0 0 ~)'• 0 o .o.o • o O0
30 20 10 0
o
1'o
2'o 3'0 4'0
5'0
6'0 7'0 8'o 90
E A S T (miles)
Fig. 6. Potential and existing sampling points (solid circles are existing data points). 100
@@ ® @ @®@® @®@® @@®@ ®@@@@ @®®®® @@®® @Q
90 80 70
60"~ "l- 50 nO 4O z 30 20 10
0
1'0
2'0 3'0 4'o
5'0
6'0 7'o
EAST (miles) Fig. 7. Variance reduction analysis for leakance.
go
9o
98
200-
(l/dr
OTV
100
)
1'0
2'0
3'0
RANK Fig. 8. TOTV and TVR at each round of kriging for leakance. higher priorities must be assigned to these zones when planning a monitoring network. The central zone and the western boundary, on the other hand, have lower ranks which means that they should not be considered as primary sampling locations. Figure 8 illustrates the level of u n c e r t a i n t y at each round of sampling in terms of T O T V (i.e., total sum of estimation variances), and the amount of information gain due to each new sampling in terms of its T V R , eqn. (7). As shown, sampling at the top five points yields the highest amount of information gain, while additional measurements appear to produce only marginal benefits. One can assume that there must be a finite number of new sampling sites, such that any additional sampling will result in information gains that cannot be justified in economic terms. S a m p l i n g based on median r a n k i n g
Median ranking, shown in Fig. 9, provides an alternative criterion based on the estimated magnitudes of leakance. However, the accuracy of the estimated
99 100
70
@@ ® @ ®@@@
60
@@@@
90 80
E "~ l- 50 nO 40 z
@@®® ®@@®@ @@@®® ®®@Q
30 20
®@
10 0
1'0
20
3'0
4'0
5'0
6'0
7'0
8'0
90
EAST (miles) Fig. 9. Median ranking for leakance.
100"
90
® ®
80
®
70
@®®@
~6o
@@@®
l.-rr
®@@@@
040 z
30
® ® @ @ ®
20
®@@Q
,o
•
@
20
3'0
0
@
0
1'0
4'0
5'0
EAST (miles) Fig. 10. 10% risk ranking for leakance.
6'0
7'0
8'0
90
1O0
100
O@ O0 @0®@ 0@00 000(~ 00000 Q 0000 00@@ QQ
90 80 70 60 E v 50 -r knO 40 Z 30 20 10 0
1'0
2'0
3'0
4'0
5'0
6'0
7~0
8'0
90
EAST (miles) Fig. 11. 1%, risk ranking for leakance.
100-
O0 0 0 O@~)Q 0@00 0@00 00000 ®000© 00@@ @@
90 80 70 60 ~
50
0 z
40
I
30 20 10
0
10
2'o 3'0 4'0
s'o 0'0 7'0 8'0 90
EAST (miles) Fig. 12.40%, risk ranking for leakance.
101
field has been ignored. In this ranking the southern tip and the upper central portions of the Dougherty Plain have the highest priorities, while the lower central and the n o r t h e r n tip zones and the eastern boundary have acquired lower priorities. Sampling based on risk ranking Risk value which depends on both the magnitude and the accuracy of an estimated leakance, seems to be a more suitable criterion for our sampling purposes. As illustrated in Fig. 10, the southern tip and the central zones have been identified as the primary candidates for sampling. The risk ranking, unlike the variance reduction analysis, has not ignored areas of high leakage with good accuracy, such as the central region (nodes ranked 4th, 5th, and 6th). On the other hand, areas with low potential recharge, such as the lower central and the n o r th er n tip zones along with the western boundary, have received lower priorities. These results demonstrate the advantages of risk ranking, which in addition to the variance, considers the actual magnitudes of leakance, and yields a more balance guideline for sampling. Equation (8) shows that the risk value is basically a weighted sum of the expected value and the estimation variance of the log-transformed leakance. Thus, as we decrease the probability of exceedence of our risk values, we are in fact giving more weight to the variance, and conversely. For instance, the 1% risk ranking, shown in Fig. 11, displays a move toward the results of variance reduction analysis as shown in Fig. 7. On the other hand, as expected, the result of the 40% risk ranking, shown in Fig. 12, is almost identical to the median ranking given in Fig. 9. The risk ranking gives us a degree of flexibility to adjust this weighted sum according to our specific sampling objectives. This work clearly demonstrates the potentials of geostatistical techniques in the design of regional groundwater monitoring networks. Following these optimally oriented procedures can lead to more efficient sampling schemes, at lower costs. ACKNOWLEDGMENTS This research was supported by grant ECE-8503897 from the National Science Foundation, and U.S. Department of Interior, Geological Survey Project G-1219(05). The authors would like to thank Prof. Ram Arora of Georgia State University for his help and guidance.
REFERENCES Benjamin, J.R. and Cornell, C.A., 1970. Probability, Statistics, and Decision for Civil Engineers. McGraw-Hill, New York, N.Y. Delfiner, P., 1975. Linear estimation of non-stationary spatial phenomena. In: M. Guarascio. M. David, and C. Huijbregts (Editors), AdvancedGeostatistics in the Mining Industry. Nato Adv. Stud. Inst.. Reidel, Hingham, Mass.
102 Freeze, R.A., 1975. A stochastic-conceptual analysis of one-dimensional groundwater flow in nonuniform homogeneous media. Water Resour. Res., 11(5): 72~741. Georgia Department of Agriculture, 1985. Georgia Agricultural Facts. Ga. Dep. Agric., Atlanta, Ga. Hayes, L.R., Maslia, M.L. and Meeks, W.C., 1983. Hydrology and model evaluation of the principal artesian aquifer, Dougherty Plain, Southwest Georgia. Bull. 97, Dep. Nat. Resour.. Environ. Prot. Div., Ga. Geol. Surv., Atlanta, Ga. Journel, A.G. and Huijbregts, C., 1978. Mining Geostatistics. Academic Press, London, 600 pp. Kaf?itsas, J. and Bras, R.L., 1981. The practice of kriging. Rep. No. 263, Ralph M. Parsons I,ab.. M.I.T., Cambridge, Mass. McConnel, J.B., Hicks, D.W., Lowe, L.E., Cohen, S.Z. and Jovanovich, A.P., 1984. Investigation of ethylene dibromide (EDB) in ground water in Seminole County, Georgia. U.S. Dep. Inter., Geol. Surv., Circ. 933, Reston, Va. Meyer, S.T., 1975. Data Analysis for Scientists and Engineers. Wiley, New York, N.Y., 285 pp. Rouhani, S.. 1983. Optimal data collection in random fields. Ph.D. Thesis, Div. Appl. Sci., Harvard Univ., Cambridge, Mass. Rouhani, S., 1985. Variance reduction analysis. Water Resour. Res., 21(6): 837-846. Rouhani, S. and Fiering, M.B., 1986. Resilience of a statistical sampling scheme. J. Hydrol., 89: 1 11. Rouhani, S. and Hall, T.J., 1987. Optimal schemes for ground water quality monitoring in the shallow aquifer, Dougherty Plain, southwestern Georgia. Tech. Complet. Rep., USDI/USGS Proj. G-1219 (05), ERC 03 87, Environ. Resour. Center, Georgia Inst. Technol., Atlanta, Ga., 49 pp.