DEVELOPMENT OF GUIDELINES FOR DESIGN OF SAMPLING PROGRAMS TO PREDICT GROUNDWATER DISCHARGE by LORIE S. C A H N B.A., The University of California at Santa Cruz, 1979 B.A., The University of California at Santa Cruz, 1979 A THESIS SUBMITTED IN P A R T I A L F U L F I L L M E N T OF T H E REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF SCIENCE in T H E F A C U L T Y OF G R A D U A T E STUDIES Department of Geological Sciences We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH C O L U M B I A February 1987 ® LORIE S. C A H N , 1987 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at The University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Geological Sciences The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date: February 1987 ABSTRACT The objective of this study is to develop guidelines for the design of sampling programs to predict groundwater discharge. A method for choosing a preferred sampling strategy from a set of alternatives is presented. A framework is outlined, in the form of an objective function, that incorporates both the cost of collecting data and the worth of data. A monetary value is assigned to the worth of hydraulic conductivity data by examining the economic losses associated with the uncertainty in predictions of groundwater discharge. The method is applied to the problem of designing a sampling program that measures hydraulic conductivity for predicting discharge from a rapid infiltration pond. Hydraulic conductivitj' data are generated for hypothetical hillslopes using a stochastic finite element model. A set of sampling strategies are selected. For each sampling strategy, the value and location of measurements and the uncertainty in the spatial variation of hydraulic conductivity are incorporated using conditional simulations. Estimates of pond discharge are calculated from the stream function solution and compared to the actual value of pond discharge for the hypothetical site. The root mean square error is used to quantify the uncertainty in discharge predictions. A set of alternative sampling strategies are compared using the objective function. Prediction uncertainty, measured by the root mean square error, is sensitive to both the structure of the heterogeneities and the location of measurements. Sampling schemes that lead to good estimates of the ensemble mean and standard deviation will not necessarily lead to good predictions of discharge. The goal of sampling schemes should be to collect data in key locations of the flow domain and to identify the spatial variation in hydraulic i i conductivity in a cost effective manner. For predicting discharge from a rapid infiltration pond, locating one or two initial boreholes below the pond is the preferred strategy for a majority of the cases tested. When the measurements are spaced evenly throughout the flow domain, important shallow layers may be missed that have a large influence on pond discharge. Increasing the number of boreholes does not necessarily lead to more certain predictions of pond discharge or to lower values of the objective function. Considerable uncertainty in discharge predictions can exist even with a relatively large number of measurements. While an optimal strategy exists, there is potential for significant variation in prediction uncertainty at individual sites. iii T A B L E O F C O N T E N T S ABSTRACT i i L i s t of Tables v i L i s t of F i g u r e s v i i ACKNOWLEDGEMENTS v i i i INTRODUCTION 1 PREVIOUS WORK • 3 S p a t i a l v a r i a b i l i t y of h y d r a u l i c c o n d u c t i v i t y .... 3 S t o c h a s t i c models 4 Network design 11 PURPOSE 17 OBJECTIVES 19 METHODOLOGY 22 STOCHASTIC MODEL 27 S p a t i a l v a r i a b i l i t y 27 Generating the re f e r e n c e case 29 C o n d i t i o n a l s i m u l a t i o n 31 SOLUTION 35 Stream f u n c t i o n s 35 Boundary c o n d i t i o n s 37 Discharge 37 EVALUATION OF SAMPLING STRATEGIES 38 Data c o l l e c t i o n c o s t s 38 Data worth 40 O b j e c t i v e f u n c t i o n 46 SOURCE OF COMPUTER CODES 46 COMPUTER SIMULATIONS . 48 PHYSICAL MODEL 48 Boundary value problem 48 F i n i t e element g r i d 49 STOCHASTIC SIMULATIONS 50 Reference cases 50 S p a t i a l averaging 52 Sample g r i d s 52 C o n d i t i o n a l s i m u l a t i o n input 60 C o n d i t i o n a l s i m u l a t i o n output 62 COST COEFFICIENTS 63 COMPUTER TIME 64 RESULTS 66 HOMOGENEOUS CASE 66 A PRIORI (UNCONDITIONAL) MODEL 68 REFERENCE CASES 1 THROUGH 5 75 iv LOCATION FOR INITIAL BOREHOLES 89 Sr used as input 92 Sy used as input 96 IMPORTANT REGIONS TO SAMPLE 100 Reference cases with small ay 101 Reference cases with l a r g e r standard d e v i a t i o n . 109 Reference cases with s m a l l e r i n t e g r a l s c a l e s ... 117 LOCATION FOR SECOND BOREHOLE 122 Reference cases 11 through 13 122 SAMPLING IN STRATIFIED MEDIA 124 MULTIPLE BOREHOLES 134 OBJECTIVE FUNCTION 137 LIMITATIONS 146 SUMMARY 149 CONCLUSIONS 149 RECOMMENDATIONS 151 BIBLIOGRAPHY 153 APPENDIX : TABLE OF NOTATION 158 v Li s t of Tables 1. Input and output f o r r e f e r e n c e cases 51 2. C o n d i t i o n a l s i m u l a t i o n s 60 3. S i m u l a t i o n s with 200, 300, and 400 r e a l i z a t i o n s 61 4. Cost c o e f f i c i e n t s 64 5. S e n s i t i v i t i e s f o r di s c h a r g e s t a t i s t i c s 69 6. S and the number of independent zones 74 7. I n f l u e n c e of Y on Q{ 77 8. S i m u l a t i o n s A-D f o r r e f e r e n c e cases 1-5 using Sr 90 9. S i m u l a t i o n s A-D f o r r e f e r e n c e cases 1 - 5 using^ s ' " ....97 10. Summary of r e s u l t s f o r i n i t i a l boreholes .....100 11. S i m u l a t i o n s E-H f o r r e f e r e n c e cases 1-5 102 12. S i m u l a t i o n s E-H f o r r e f e r e n c e cases 6-10 115 13. S i m u l a t i o n s E-H f o r r e f e r e n c e cases 11-13 121 14. S i m u l a t i o n s I-K f o r r e f e r e n c e cases 11-13 123 15. S i m u l a t i o n s M-0 f o r r e f e r e n c e cases 14-17 131 16. S i m u l a t i o n s P-R f o r r e f e r e n c e cases 6-10 135 17. O b j e c t i v e f u n c t i o n f o r the s i m u l a t i o n s 139 v i Li st of Fi gures 1. V e r t i c a l s e c t i o n i l l u s t r a t i n g c e n t r a l concepts 18 2. Water t a b l e p o s i t i o n s 25 3. S t o c h a s t i c b l o c k s 30 4. F u n c t i o n e x p r e s s i n g the square of the b i a s 41 5. Flow domain and boundary c o n d i t i o n s 49 6. F i n i t e element g r i d 50 7. I n i t i a l boreholes A through D 54 8. L o c a t i o n of measurements i n regions E - H 55 9. Sampling schemes I - L with two boreholes 56 10. Sampling s t r a t e g i e s M - 0 57 11. Sampling s t r a t e g i e s P - R 59 12. Streamline and f l u x p l o t s f o r the homogeneous case ...67 13. Frequency d i s t r i b u t i o n s f o r the a p r i o r i model 71 14. Discharge from the a p r i o r i model and cases 1 - 5 ....78 15. Reference case 1 80 16. Reference cases 2 - 5 86 17. Discharge histograms from sampling regions E - H f o r case 1 104 18. Reduction of u n c e r t a i n t y i n K due to measurements ...108 19. Log K maps and s t r e a m l i n e p l o t s f o r cases 6 - 10 ....110 20. Flux p r o f i l e s f o r r e f e r e n c e cases 6 - 10 111 21. Discharge f o r the a p r i o r i model and cases 6 - 10 ...113 22. Reference cases 11 - 13 118 23. Discharge from the a p r i o r i model and cases 11-13 ...120 24. Reference cases 14 - 17 126 25. Discharge from the a p r i o r i model and cases 14-17 ...128 26. Streamline p l o t s f o r r e f e r e n c e cases 15 and 16 129 v i i ACKNOWLEDGEMENTS This research was supported by a strategic grant from the Natural Sciences and Engineering Research Council of Canada. I wish to thank my committee members, A l Freeze and Bill Caselton, and especially Les Smith, who conceived this project and with whom I worked closely. I feel very fortunate to have been part of the Freeze-Smith connection. Sincere thanks to Craig Forster for numerous brainstorming sessions, Joel Massmann and Reidar Zapf'-Gilje for useful discussions, Chuck Mase, J im Glosli, and Tom Nicol for sharing their useful tricks for cajoling the Array Processor, Gord Hodge for drafting the figures, and Larry Roberts for helping to TeX the tables. A personal note of thanks to my father, John, who gave up precious time during visits with his grandson to discuss my research, and to Doug who patiently survived the West Coast rain just a little longer than anticipated. viii INTRODUCTION S t o c h a s t i c methods pr o v i d e a p r o b a b i l i s t i c framework fo r d e s c r i b i n g the n a t u r a l v a r i a b i l i t y of g e o l o g i c media. In hydrogeology, these methods have been p r i m a r i l y used to i n v e s t i g a t e the p h y s i c s of flow and t r a n s p o r t i n heterogeneous media or to q u a n t i f y u n c e r t a i n t y i n model p r e d i c t i o n . The use of these models in developing g u i d e l i n e s for the design of sampling programs has r e c e i v e d l i t t l e a t t e n t i o n . The goal of t h i s study i s to determine whether key l o c a t i o n s e x i s t f o r c o l l e c t i n g data that tend to reduce the u n c e r t a i n t y in model p r e d i c t i o n . In a d d i t i o n , a simple framework i s intr o d u c e d f o r e v a l u a t i n g the worth of data c o l l e c t e d . Using t h i s framework, the t r a d e o f f between the cost of c o l l e c t i n g samples and the r e d u c t i o n i n p r e d i c t i o n u n c e r t a i n t y can be q u a n t i f i e d . F i n a l l y , t h i s framework i s used to i n v e s t i g a t e whether q u a l i t a t i v e g u i d e l i n e s f o r the design of sampling networks can be determined. T h i s i s a conceptual study, s i m i l a r to what Freeze [1980] r e f e r r e d to as a " h y d r o l o g i c game". I t i s a game i n the sense that we work with h y p o t h e t i c a l h i l l s l o p e s f o r which we have p e r f e c t knowledge of the u n d e r l y i n g hydrogeologic parameters and pro c e s s e s . We p l a y the game by askin g "what i f " - t y p e q u e s t i o n s . What i f we gathered only so much in f o r m a t i o n upon which to base our d e c i s i o n s ? How wrong would our p r e d i c t i o n s have been and what would be the 1 INTRODUCTION / 2 economic consequences of being wrong? It i s hoped that t h i s study w i l l a i d in d e v e l o p i n g sampling s t r a t e g i e s f o r p r e d i c t i n g groundwater d i s c h a r g e . T h i s study i s not, however, an o p e r a t i o n a l t o o l f o r s p e c i f i c f i e l d s i t e s . The methodology developed here can not be used at an a c t u a l f i e l d s i t e t o design an optimum sampling program. Instead, the aim of t h i s study i s to i n v e s t i g a t e g u i d e l i n e s f o r the design of sampling programs by examining the e f f e c t i v e n e s s of a f i n i t e number of sampling s t r a t e g i e s in reducing p r e d i c t i o n u n c e r t a i n t y at d i f f e r e n t h y p o t h e t i c a l s i t e s . U n c e r t a i n t y i n p r e d i c t i o n s from groundwater models stem from s e v e r a l sources of u n c e r t a i n t y in input parameters. Hydrogeologic data c o n t a i n two types of u n c e r t a i n t y : i n t r i n s i c and i n f o r m a t i o n u n c e r t a i n t y [Detlinger and Wilson, 1981]. I n t r i n s i c u n c e r t a i n t y stems from the s p a t i a l v a r i a b i l i t y of hydrogeologic parameters due to n a t u r a l v a r i a t i o n s i n g e o l o g i c p r o c e s s e s . T h i s i n h e r e n t type of u n c e r t a i n t y occurs at a l l s c a l e s and i s i r r e d u c i b l e . Even with a l a r g e number of measurements, the true h e t e r o g e n e i t y can not be p r e c i s e l y c h a r a c t e r i z e d . A groundwater model i s , n e c e s s a r i l y , only a " . . . s i m p l i f i e d v e r s i o n of complex r e a l i t y " [Smith and Schwartz, 1981 a ]. Rather than attempting to i d e n t i f y the a c t u a l s p a t i a l arrangement of parameters, t h i s v a r i a b i l i t y can be d e s c r i b e d i n a p r o b a b i l i s t i c framework through the use of s t o c h a s t i c models. INTRODUCTION / 3 Information u n c e r t a i n t y stems from e r r o r s due to noisy or i n s u f f i c i e n t data and i s r e d u c i b l e through e f f e c t i v e sampling s t r a t e g i e s [Det l i n g e r and Wilson, 1981]. T h i s study focuses on the i d e n t i f i c a t i o n of key l o c a t i o n s i n a flow f i e l d where measurements can reduce i n f o r m a t i o n u n c e r t a i n t y , thereby reducing p r e d i c t i o n u n c e r t a i n t y . PREVIOUS WORK S p a t i a l v a r i a b i l i t y of h y d r a u l i c c o n d u c t i v i t y The n a t u r a l v a r i a b i l i t y of a q u i f e r m a t e r i a l s c r e a t e s p a t t e r n s of groundwater flow that may d i f f e r s i g n i f i c a n t l y from those in uniform porous media. The accuracy of a p r e d i c t i o n from a groundwater model i s a f f e c t e d by our a b i l i t y to map the v a r i a t i o n inherent i n a q u i f e r p r o p e r t i e s , such as h y d r a u l i c c o n d u c t i v i t y . Due to the complexity of n a t u r a l v a r i a t i o n at a l l s c a l e s , i t i s g e n e r a l l y impossible to i d e n t i f y the a c t u a l v a r i a t i o n s f o r use i n a d e t a i l e d d e t e r m i n i s t i c model. The t r e n d i n the l a s t decade has been to d e s c r i b e these v a r i a t i o n s using a p r o b a b i l i s t i c framework. S e v e r a l authors r e p o r t the r e s u l t s of d e t a i l e d measurements of h y d r a u l i c c o n d u c t i v i t y at v a r i o u s s i t e s to c h a r a c t e r i z e the s t a t i s t i c a l p r o p e r t i e s [ F r e e z e , 1975; Smith, 1981; El Kadi and Brutsaert, 1985]. For a comprehensive l i s t of f i e l d s t u d i e s a d d r e s s i n g the s p a t i a l INTRODUCTION / 4 v a r i a b i l i t y of h y d r o l o g i c parameters, the reader i s r e f e r r e d to Loague [1986]. S t o c h a s t i c models The o b j e c t i v e of s t o c h a s t i c models i s to i n c o r p o r a t e the u n c e r t a i n t y i n input parameters by d e s c r i b i n g them s t a t i s t i c a l l y r a t her than d e t e r m i n i s t i c a l l y . The f i r s t s t o c h a s t i c model of flow through heterogeneous media was developed by Warren and P r i c e [1961]. Values of h y d r a u l i c c o n d u c t i v i t y , chosen randomly from a frequency d i s t r i b u t i o n , were as s i g n e d to d i s c r e t e b l o c k s used i n the numerical s o l u t i o n of the flow e q u a t i o n . Using a Monte C a r l o method, the process was repeated with new values of h y d r a u l i c c o n d u c t i v i t y u n t i l a l a r g e number of s o l u t i o n s were obtained that c o u l d be analysed s t a t i s t i c a l l y . T h i s p i o n e e r i n g work i n v e s t i g a t e d the i n f l u e n c e of d i f f e r e n t frequency d i s t r i b u t i o n s f o r h y d r a u l i c c o n d u c t i v i t y , i n c l u d i n g lognormal, on steady and t r a n s i e n t flow. The work of V/arren and P r i c e was extended by Freeze [1975] to i n c l u d e p o r o s i t y , c o m p r e s s i b i l i t y , and h y d r a u l i c c o n d u c t i v i t y as input parameters f o r a one-dimensional model. The r e s u l t s l e d Freeze to q u e s t i o n the v a l i d i t y of r e p r e s e n t i n g heterogeneous media with a homogeneous model. Se v e r a l of h i s assumptions, however, were u n r e a l i s t i c . F l u i d s i n one-dimensional problems cannot flow around blocks INTRODUCTION / 5 of low h y d r a u l i c c o n d u c t i v i t y . A d d i t i o n a l l y , although the frequency d i s t r i b u t i o n s of the three input parameters were c r o s s - c o r r e l a t e d , the values a s s i g n e d to each parameter were s t a t i s t i c a l l y independent. One would expect h y d r a u l i c p r o p e r t i e s to d i s p l a y some degree of s p a t i a l p e r s i s t e n c e . The work was extended by Smith and Freeze [1979a] to i n c l u d e s p a t i a l a u t o c o r r e l a t i o n u sing a nearest neighbor model and was l a t e r extended to two dimensions [Smith and Freeze, 19796]. A drawback of the nearest neighbor method i s that i t does not handle c o v a r i a n c e f u n c t i o n s e x p l i c i t l y . S e v e r a l techniques have been developed f o r g e n e r a t i n g random f i e l d s where the s p a t i a l s t r u c t u r e can be s p e c i f i e d as i n p u t . The technique based on s p e c t r a l a n a l y s i s generates continuous random f i e l d s . Using s p e c t r a l methods, Mejia and Rodri guez-It urbe [1974] developed a procedure f o r s y n t h e s i z i n g i s o t r o p i c random f i e l d s . A more e f f i c i e n t technique f o r g e n e r a t i n g s p a t i a l v a r i a b i l i t y i n l a r g e domains, the t u r n i n g bands method, was developed by Mat heron [1973], The t u r n i n g bands method uses s p e c t r a l techniques to transform m u l t i d i m e n s i o n a l problems i n t o a sum of a s e r i e s of one-dimensional l i n e p r o c e s s e s . O r i g i n a l l y , i t was a p p l i e d to three-dimensional mining problems [ J o u r n e l and Hui j br egt s, 1978] and l a t e r to the two-dimensional case. The two-dimensional case i s more d i f f i c u l t to simulate than the three-dimensional case because i t s one-dimensional INTRODUCTION / 6 eq u i v a l e n c e i s more complicated [Mantoglou and Wilson, 1982], Mantoglou and Wilson g e n e r a l i z e d t h i s approach to i n c l u d e any covariance f u n c t i o n [1982] and extended i t to a n i s o t r o p i c and a r e a l averaged processes [1981]. A much simpler method of g e n e r a t i n g d i s c r e t e values of a s p a t i a l l y c o r r e l a t e d random f i e l d [Scheuer and Stol I er , 1962; Clifton and Neuman, 1982] i s used i n t h i s study and w i l l be d i s c u s s e d l a t e r . Random f i e l d s can be used as input to numerical groundwater models by s u b d i v i d i n g the flow domain i n t o d i s c r e t e regions and a s s i g n i n g a s i n g l e value of the random f i e l d to each r e g i o n . In the Monte C a r l o approach, many r e a l i z a t i o n s of a s t o c h a s t i c process are generated and the p r o b a b i l i t y d i s t r i b u t i o n s of the output v a r i a b l e s c a l c u l a t e d . The advantage of t h i s technique i s that i t can be a p p l i e d to boundary value problems with complicated geometries, boundary c o n d i t i o n s and input parameters. The main disadvantages of t h i s technique are that i t r e q u i r e s l a r g e core storage and e x c e s s i v e amounts of computing time. In order to preserve the s p a t i a l s t r u c t u r e of the random f i e l d s , the d i s c r e t i z e d regions must be small r e l a t i v e to the average d i s t a n c e over which a parameter i s c o r r e l a t e d . T h i s d i s t a n c e , r e f e r r e d to as the i n t e g r a l s c a l e , must be small r e l a t i v e to the l e n g t h of the flow domain being modeled. For these reasons the Monte C a r l o approach r e q u i r e s INTRODUCTION / 7 dense g r i d s and l a r g e m a t r i c e s . With the advances i n computer technology of the l a s t decade, these l i m i t a t i o n s are becoming l e s s of a concern. There are s e v e r a l a l t e r n a t i v e s to the Monte C a r l o method f o r s o l v i n g s t o c h a s t i c flow equations. In p e r t u r b a t i o n methods, only the f i r s t two moments (mean and covar i a n c e ) of the output v a r i a b l e s are d e r i v e d . Using s p e c t r a l a n a l y s i s , the input parameters can be w r i t t e n i n terms of a mean and a p e r t u r b a t i o n . An equation i s d e r i v e d and s o l v e d a n a l y t i c a l l y that l i n k s the p e r t u r b a t i o n i n h y d r a u l i c c o n d u c t i v i t y (K) with the p e r t u r b a t i o n i n h y d r a u l i c head (h) [Bakr et al., 1978]. The main advantages of t h i s method are that i t rep r e s e n t s the s p a t i a l s t r u c t u r e of the porous medium as a continuum and that i t pr o v i d e s an economical, c l o s e d form s o l u t i o n f o r the s t a t i s t i c a l p r o p e r t i e s of h y d r a u l i c head. The major disadvantage of s p e c t r a l s o l u t i o n s i s that the output d i s t r i b u t i o n s must be s t a t i o n a r y . T h i s r e q u i r e s that the s t a t i s t i c a l moments be constant i n space. For t h i s reason, the s p e c t r a l method i s l i m i t e d to very simple geometries, boundary c o n d i t i o n s and input d i s t r i b u t i o n s . I t i s not d i r e c t l y a p p l i c a b l e to f i n i t e domains. In a d d i t i o n , the v a r i a n c e of the input parameter must be s m a l l . L a s t l y , the p r e s e r v a t i o n of f i e l d measurements i s not s t r a i g h t f o r w a r d [Dagan, 1982a]. A numerical p e r t u r b a t i o n method presented by Sagar INTRODUCTION / 8 [1978] i n v o l v e s T a y l o r s e r i e s expansions of the s o l u t i o n of the governing equation around expected val u e s of K and h [ D e t l i n g e r and Wilson, 1981; Townley, 1984], T h i s approach i s capable of d e a l i n g with nonuniform flow, t r a n s i e n t c o n d i t i o n s , and complex geometries and boundary c o n d i t i o n s . The main disadvantage i s that i t i s r e s t r i c t e d to input v a r i a n c e s that are s u f f i c i e n t l y small to ensure that the c o e f f i c i e n t of v a r i a t i o n i s a f r a c t i o n of one [Det l i n g e r and Wilson, 1981]. The techniques f o r generating random f i e l d s d e s c r i b e d above are c a l l e d u n c o n d i t i o n a l because the l o c a t i o n of data p o i n t s i s not c o n s i d e r e d . U n c o n d i t i o n a l models have been used p r i m a r i l y f o r i n v e s t i g a t i n g the p h y s i c s of flow i n heterogeneous media. In these models, the v a r i a n c e of the input parameters, a measure of the u n c e r t a i n t y i n the estimate at any l o c a t i o n , i s constant i n space. If f i e l d measurements of input parameters are a v a i l a b l e , the u n c e r t a i n t y i n these parameters can be reduced i n the v i c i n i t y of the measurements. The technique f o r p r e s e r v i n g f i e l d data i s c a l l e d c o n d i t i o n a l s i m u l a t i o n and was developed by Mat heron [1973]. In c o n d i t i o n a l s i m u l a t i o n s , a l a r g e number of r e a l i z a t i o n s of the random f i e l d are generated in which the measurements are preserved, but the remainder of the f i e l d f l u c t u a t e s randomly a c c o r d i n g to s p e c i f i e d s t a t i s t i c s INTRODUCTION / 9 d e s c r i b i n g the s p a t i a l v a r i a b i l i t y . In t h i s way random f i e l d s are generated which are p o s s i b l e v e r s i o n s of the true f i e l d because they are c o n s i s t e n t with the measurements. Del homme [1979], the f i r s t h y d r o l o g i s t to use t h i s technique, mapped h y d r a u l i c heads in the Bathonian a q u i f e r in France. He found that p r e s e r v i n g measurements of t r a n s m i s s i v i t y in the c o n d i t i o n a l model was not always s u f f i c i e n t to ensure a good f i t with h y d r a u l i c head data and only m a r g i n a l l y reduced the u n c e r t a i n t y (variance) i n p r e d i c t e d heads. The e f f e c t of measurements depended upon the l o c a t i o n of the measurements w i t h i n the flow regime, e s p e c i a l l y with re s p e c t to the boundaries of the system. Delhomme suggested that measuring t r a n s m i s s i v i t y in key l o c a t i o n s i n the flow regime would improve the q u a l i t y of p r e d i c t i o n s . T h i s i m p l i e s that the u s e f u l n e s s of a measurement for making p r e d i c t i o n s depends upon the l o c a t i o n . Key l o c a t i o n s are those l o c a t i o n s where the value of i n f o r m a t i o n gained from, a measurement i s h i g h . C l i f t o n and Neuman [1982] demonstrated that the e f f e c t of c o n d i t i o n i n g c o u l d be q u i t e important i n reducing the u n c e r t a i n t y in p r e d i c t i o n s of h y d r a u l i c head. T h e i r data, i n comparison with Delhomme's, were more evenly and densely d i s t r i b u t e d throughout the flow domain. In a d d i t i o n , the boundary value problems d i f f e r e d , p a r t i c u l a r l y with respect to boundary c o n d i t i o n s . C l i f t o n and Neuman f u r t h e r improved INTRODUCTION / 10 t h e i r p r e d i c t i o n s by using an i n v e r s e model to i n c o r p o r a t e measurements of both h y d r a u l i c head and flow r a t e i n a d d i t i o n to the t r a n s m i s s i v i t y d a t a . C o n d i t i o n a l s i m u l a t i o n s have a l s o been used to i n v e s t i g a t e the r o l e of h y d r a u l i c c o n d u c t i v i t y data in reducing the u n c e r t a i n t y i n p r e d i c t i o n s of mass t r a n s p o r t . Smith and Schwartz [1981 a] found that the l o c a t i o n of zones of higher h y d r a u l i c c o n d u c t i v i t y dominated the p a t t e r n s of mass t r a n s p o r t and suggested that i d e n t i f y i n g these zones should be e s s e n t i a l f o r a c c u r a t e p r e d i c t i o n s . T h e i r subsequent work [19816] focused on how the number of data p o i n t s a f f e c t s p r e d i c t i o n u n c e r t a i n t y . They demonstrated that even q u i t e a l a r g e number of h y d r a u l i c c o n d u c t i v i t y measurements do not go f a r toward reducing the u n c e r t a i n t y i n t r a n s p o r t p r e d i c t i o n s . U n c e r t a i n t i e s in p r e d i c t i o n were r e l a t e d i n a complex way to such f a c t o r s as the number of data p o i n t s , the c o r r e l a t i o n d i s t a n c e , the l o c a t i o n of data p o i n t s with respect to zones of higher h y d r a u l i c g r a d i e n t and the standard d e v i a t i o n i n h y d r a u l i c c o n d u c t i v i t y . They conclude that "... c o n s i d e r a b l e h y d r a u l i c c o n d u c t i v i t y data may be necessary to o b t a i n a reasonable degree of confidence i n p r e d i c t i o n s of s i t e behavior". INTRODUCTION / 11 Network design The above s t u d i e s use s t o c h a s t i c models p r i m a r i l y to i n v e s t i g a t e the p h y s i c s of flow and t r a n s p o r t i n heterogeneous media, and to q u a n t i f y the u n c e r t a i n t y in model p r e d i c t i o n . A l o g i c a l next step i s to use t h i s i n f o r m a t i o n i n d e s i g n i n g data c o l l e c t i o n networks. S e v e r a l r e s e a r c h e r s have touched upon t h i s t o p i c . T h e i r approaches depend upon the s p e c i f i c o b j e c t i v e s of the sampling program. The o b j e c t i v e s focus on minimizing one or more of the f o l l o w i n g sources of p r e d i c t i o n u n c e r t a i n t y or e r r o r (modified from [Rodriguez-Iturbe and Mejia, 1974; Bogardi, et a l . , 1985]) Level 1 e r r o r s i n p r e d i c t i n g input parameters, such as h y d r a u l i c c o n d u c t i v i t y ; Level 2 e r r o r s i n p r e d i c t i n g output, such as h y d r a u l i c head or d i s c h a r g e , r e l a t e d to l e v e l 1 through p h y s i c a l processes such as groundwater flow; Level 3 net economic l o s s r e l a t e d to d e c i s i o n s made based on l e v e l 1 or l e v e l 2 e r r o r s . M i n i m i z i n g l e v e l 1 e r r o r s i s a problem of s t a t i s t i c a l i n f e r e n c e [Dagan, 1985]. Data worth i s i n t e r p r e t e d i n terms of the r e d u c t i o n i n u n c e r t a i n t y i n e s t i m a t i n g a parameter. INTRODUCTION / 12 The c o s t of c o l l e c t i n g data may be c o n s i d e r e d but the economic value of the data i s not. In g e n e r a l , two approaches have been f o l l o w e d . A network can be designed to minimize the e r r o r of e s t i m a t i o n w i t h i n a f i x e d budget [Hughes and Lettenmaier, 1981]. An a l t e r n a t i v e i s to minimize sampling c o s t s s u b j e c t to the c r i t e r i o n of minimal a c c e p t a b l e accuracy. These c r i t e r i a are computed using v a r i a n c e r e d u c t i o n f a c t o r s or k r i g i n g . K r i g i n g i s used to c a l c u l a t e the e s t i m a t i o n v a r i a n c e , a measure of the u n c e r t a i n t y of the estimate, f o r p o t e n t i a l measurement p o i n t s . The e s t i m a t i o n v a r i a n c e can be computed f o r a p o t e n t i a l p o i n t without a c t u a l l y t a k i n g a measurement at that p o i n t . The only i n f o r m a t i o n necessary i s the c o v a r i a n c e f u n c t i o n or variogram and the geometry of the measurement p o i n t s . The o b j e c t i s to sample in such a way that the e s t i m a t i o n v a r i a n c e i s minimized. K r i g i n g [Delhomme, 1983] and v a r i a n c e r e d u c t i o n f a c t o r s [Rodriguez-It urbe and Mejia, 1974] have been used to optimize the l o c a t i o n of an a d d i t i o n a l measurement p o i n t . T h i s s e q u e n t i a l method pr o v i d e s an optimal s o l u t i o n f o r an a d d i t i o n a l p o i n t but i s not n e c e s s a r i l y optimal f o r a set of measurements. In order to be the optimal design f o r a set of p o i n t s , the p o i n t s would have to be s e l e c t e d s i m u l t a n e o u s l y . A technique has been proposed f o r augmenting an e x i s t i n g network with one or more p o i n t s based on a "best" design INTRODUCTION / 13 i n s t e a d of an optimum design [Bogardi et al . , 1985]. In t h i s sense, a best design i s one which i s optimal f o r a subset of d e s i g n s . B e t t e r designs may e x i s t , such as the optimum des i g n , but they are not members of the subset of designs c o n s i d e r e d . Carrera et al. [1984] combine k r i g i n g with n o n l i n e a r programming to s e l e c t an optimal set of measurement l o c a t i o n s from a f i n i t e number of p o s s i b l e p o i n t s f o r e s t i m a t i n g chemical composition. The above methods f o r network design only address l e v e l 1 e r r o r s . A groundwater flow equation i s never s o l v e d and the c o n s t r a i n t s imposed by boundary c o n d i t i o n s , or the p o s i t i o n of measurement p o i n t s i n a flow domain, are not c o n s i d e r e d . Using k r i g i n g , e s t i m a t i o n e r r o r s are c a l c u l a t e d f o r p o s s i b l e measurement p o i n t s without a c t u a l l y c o l l e c t i n g data at those p o i n t s . Combining l e v e l 1 and l e v e l 2 e r r o r s i s more complicated. I t r e q u i r e s c o l l e c t i n g data at a p o t e n t i a l measurement p o i n t f o r use i n a model of groundwater flow, and t h i s makes i t d i f f i c u l t to p r e d i c t the r e d u c t i o n i n u n c e r t a i n t y i n model output p r i o r to t a k i n g a measurement at that p o i n t . The recent r e s e a r c h combining i n v e r s e modeling with g e o s t a t i s t i c s [Neuman, 1982; Kitanidis and Vomvoris, 1983] may l e a d to network designs that combine l e v e l 1 and 2. The e f f e c t of c o l l e c t i n g data on reducing l e v e l 2 e r r o r s can be i n v e s t i g a t e d using h y p o t h e t i c a l data s e t s . The INTRODUCTION / 14 use of h y p o t h e t i c a l data s e t s , as i n t h i s study, i s u s e f u l i n i n v e s t i g a t i n g c o n c e p t u a l g u i d e l i n e s f o r network d e s i g n . Measurements from the h y p o t h e t i c a l system can be used as input to c o n d i t i o n a l models to analyse the e f f e c t of data on reducing u n c e r t a i n t i e s i n the output. Andersson et al . [1984] use a c o n d i t i o n a l model to eva l u a t e the extent to which core samples c o l l e c t e d from a h y p o t h e t i c a l f r a c t u r e d rock reduce u n c e r t a i n t y i n the geometry of the f r a c t u r e network. They suggest that at some p o i n t the b e n e f i t of making a d d i t i o n a l measurements i s low compared to the c o n s i d e r a b l e c o s t of d r i l l i n g more c o r e s . The t r a d e o f f between u n c e r t a i n t y r e d u c t i o n and the co s t of cores i s l e f t u n q u a n t i f i e d because economic value ( l e v e l 3) i s not assigned to u n c e r t a i n t y r e d u c t i o n . In many of the l e v e l 1 s t u d i e s c i t e d above, the o b j e c t i v e of the network i s to c o l l e c t base l i n e data. The economic value of base l i n e data may be d i f f i c u l t to c a l c u l a t e i n some i n s t a n c e s , which can make i t hard to combine l e v e l 1 with l e v e l 3 e r r o r s . For some types of network design, the d e s i r e to a c c u r a t e l y estimate a parameter can be l i n k e d to economic g a i n . C o l l e c t i n g data to improve an estimate i s only worthwhile i f the marginal b e n e f i t s produced by a d d i t i o n a l data are g r e a t e r than the cos t of c o l l e c t i n g the data. The t r a d e o f f between the investment i n i n f o r m a t i o n and the r e d u c t i o n i n u n c e r t a i n t y INTRODUCTION / 15 may show a d i m i n i s h i n g r e t u r n at some po i n t [Andersson et al. , 1984]. Rouhani [1985] uses 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 to s e q u e n t i a l l y s e l e c t s i t e s f o r c o l l e c t i n g measurements of h y d r a u l i c head. His approach l i n k s economic r i s k s to the u n c e r t a i n t i e s i n h y d r a u l i c head measurements by a s s i g n i n g a monetary l o s s to the l e v e l of accuracy of water l e v e l e s t i m a t e s . The amount of i n f o r m a t i o n gain i s measured by v a r i a n c e r e d u c t i o n and the expected economic gain i s measured by l o s s r e d u c t i o n . T h i s p r o v i d e s a framework f o r determining when the net expected b e n e f i t of f u r t h e r measurements becomes negative, i . e . when the c o s t of a d d i t i o n a l data exceeds the b e n e f i t of i n f o r m a t i o n g a i n . In a d d i t i o n , the e f f e c t of a new measurement on the l e v e l of accuracy of the estimated f i e l d as a whole can be c a l c u l a t e d . As a r e s u l t , the best sampling l o c a t i o n i s s e l e c t e d r e g a r d l e s s of whether the p o i n t has the h i g h e s t l e v e l of e s t i m a t i o n u n c e r t a i n t y . Because of the s e q u e n t i a l ranking of s i t e s , however, the set of s i t e s chosen i s not n e c e s s a r i l y optimal as the optimal s i t e s must be chosen si m u l t a n e o u s l y . Rouhani 's approach combines l e v e l 1 and l e v e l 3 e r r o r s , but does not address l e v e l 2 e r r o r s . H y d r a u l i c head i s regarded as an independent v a r i a b l e , r a t h e r than a f u n c t i o n of the h y d r o g e o l o g i c a l parameters. A groundwater flow INTRODUCTION / 16 equation i s never s o l v e d . Because Rouhani does not c o n s i d e r the e f f e c t of groundwater flow and boundary c o n d i t i o n s on h y d r a u l i c head, areas with low sampling d e n s i t i e s are given a c l e a r p r i o r i t y f o r sampling. In a d d i t i o n , p o i n t s on the boundary have a higher p r i o r i t y f o r sampling than i n t e r i o r p o i n t s . T h i s i s because the boundary p o i n t s are e x t r a p o l a t e d and are t h e r e f o r e l e s s r e l i a b l e than the i n t e r i o r p o i n t s , which are i n t e r p o l a t e d . In an e a r l y attempt to combine a l l three l e v e l s of e r r o r , Maddock [1973] proposed a r i s k f u n c t i o n f o r a h y p o t h e t i c a l farm to i n t e r p r e t the accuracy of data i n economic terms. He equates worth of data to the r e g r e t i n c u r r e d by l o s s i n farm income when an u n d e s i r a b l e d e c i s i o n i s made. A simple groundwater model i s combined with a management model and d e c i s i o n theory to rank d i f f e r e n t types of data. The value of the r i s k i s found to be i n s e n s i t i v e to changing the value of the h y d r o l o g i c a l parameters ( t r a n s m i s s i v i t y and storage c o e f f i c i e n t ) . Because estimates of h y d r a u l i c head from groundwater models are f a i r l y robust to changes i n t r a n s m i s s i v i t y and storage c o e f f i c i e n t , i t i s not s u r p r i s i n g that l o s s f u n c t i o n s dependent on h y d r a u l i c head are i n s e n s i t i v e t o these parameters. The recent work of Massmann and Freeze [1987a, 19876] combines Bayesian d e c i s i o n theory, r i s k - c o s t - b e n e f i t a n a l y s i s , g e o s t a t i s t i c s and c o n d i t i o n a l groundwater models INTRODUCTION / 17 to assess a l t e r n a t i v e design s t r a t e g i e s f o r waste-management f a c i l i t i e s . Data worth i s i n t e r p r e t e d i n terms of the r e d u c t i o n i n r i s k s and c o s t s and the i n c r e a s e i n b e n e f i t s to the owner-operator. A framework i s proposed to examine the t r a d e o f f s between e n g i n e e r i n g design, s i t e e x p l o r a t i o n , and m o n i t o r i n g . They do not e s t a b l i s h g u i d e l i n e s f o r s i t e e x p l o r a t i o n . They suggest that to f u l l y q u a n t i f y the value of s i t e e x p l o r a t i o n r e q u i r e s an expected-regret a n a l y s i s . PURPOSE These s t u d i e s set the stage f o r the use of s t o c h a s t i c models i n d e s i g n i n g sampling programs that combine a l l three l e v e l s of p r e d i c t i o n e r r o r . Can an economical sampling program be designed to a c c u r a t e l y p r e d i c t a dependent v a r i a b l e of monetary value? F i g u r e 1 i s a v e r t i c a l c r o s s s e c t i o n that i l l u s t r a t e s the c e n t r a l concepts i n t h i s study. A sampling program i s to be designed to c o l l e c t h y d r a u l i c c o n d u c t i v i t y data to p r e d i c t groundwater flow, Q{, through a r a p i d i n f i l t r a t i o n pond. The throughput of the pond has economic value and the accuracy of estimates of Qt are l i n k e d to economic l o s s . A piezometer nest e x i s t s i n borehole A and we wish to add another borehole at e i t h e r l o c a t i o n B or C. L e v e l 1 u n c e r t a i n t y i s h i g h e s t i n l o c a t i o n B but measurements from l o c a t i o n C used in a c o n d i t i o n a l groundwater flow model may l e a d to a more accurate INTRODUCTION / 18 pond F i g . 1. V e r t i c a l c r o s s s e c t i o n i l l u s t r a t i n g c e n t r a l concepts i n t h i s study. p r e d i c t i o n of flow through the pond ( l e v e l 2 u n c e r t a i n t y ) and a decrease i n our expected l o s s ( l e v e l 3 u n c e r t a i n t y ) . Is the monetary value of our in f o r m a t i o n gain g r e a t e r than the cost of o b t a i n i n g the a d d i t i o n a l i n f o r m a t i o n ? A framework i s developed i n t h i s study that w i l l a i d in examining these i s s u e s and in choosing a best design from a set of a l t e r n a t i v e sampling networks. The framework int r o d u c e d here i s c o n c e p t u a l . A key element of t h i s approach i s that the a c t u a l flow through the pond, Qt , i s assumed known. For the h y p o t h e t i c a l s i t e s used in t h i s study, Qt i s e a s i l y c a l c u l a t e d . At an a c t u a l f i e l d s i t e , Qt i s unknown durin g the pre-pond e x p l o r a t i o n phase. Thus, t h i s technique i s not d i r e c t l y a p p l i c a b l e to f i e l d problems. R e s u l t s from t h i s study, however, can be a p p l i e d in a q u a l i t a t i v e sense. G u i d e l i n e s f o r the design of INTRODUCTION / 19 sampling programs and the worth of hydrogeologic data are examined. In t h i s study, the data are assumed to be a b s o l u t e l y r e l i a b l e . S u b j e c t i v e hydrogeologic i n f o r m a t i o n , such as estimates of measurement e r r o r or the c o n t i n u i t y of l a y e r s between boreholes, i s not i n c o r p o r a t e d i n the a n a l y s i s . In the r e a l world, s i t e i n v e s t i g a t o r s would most l i k e l y use one of two approaches. In one approach, water l e v e l s would be recorded d u r i n g a pump t e s t , and an e f f e c t i v e h y d r a u l i c c o n d u c t i v i t y estimated f o r use i n model c a l i b r a t i o n . An a l t e r n a t i v e approach would be to i n s t a l l piezometers and measure both h y d r a u l i c c o n d u c t i v i t y and h y d r a u l i c head. H y d r a u l i c c o n d u c t i v i t y values would be a d j u s t e d i n a steady s t a t e groundwater flow model u n t i l the observed heads c o u l d be reproduced s a t i s f a c t o r i l y . T h i s r e s e a r c h d i f f e r s from p r e v i o u s work because i t i n v e s t i g a t e s how the u n c e r t a i n t y i n p r e d i c t i o n s of volume d i s c h a r g e i s a f f e c t e d by the l o c a t i o n of measurements. In a d d i t i o n , the problem of network design i s approached using a l l three l e v e l s of p r e d i c t i o n u n c e r t a i n t y . OBJECTIVES The o b j e c t i v e s of t h i s r e s e a r c h are twofold; f i r s t , to develop a framework f o r e v a l u a t i n g p o t e n t i a l sampling s t r a t e g i e s and second, to use t h i s framework to i n v e s t i g a t e INTRODUCTION / 20 g u i d e l i n e s f o r design of sampling programs. S p e c i f i c a l l y , t h i s study w i l l : 1. Develop a framework f o r examining the t r a d e o f f s between the l o c a t i o n and number of measurements, the r e d u c t i o n in u n c e r t a i n t y in output p r e d i c t i o n , the economic value of improvement in p r e d i c t i o n c e r t a i n t y , and the cost of measurements. 2. Determine r e a l i s t i c sampling c o s t s . 3. Examine how the design of a sampling network i n f l u e n c e s the u n c e r t a i n t y i n model p r e d i c t i o n . 4. Determine the u n c e r t a i n t i e s i n dependent v a r i a b l e s (output) and a method of a s s i g n i n g economic value to the u n c e r t a i n t y . 5. Define an o b j e c t i v e f u n c t i o n that can be used to i d e n t i f y the worth of a sampling network by c a l c u l a t i n g the sampling cost and the economic l o s s a s s o c i a t e d with p r e d i c t i o n u n c e r t a i n t y . 6. Determine whether key l o c a t i o n s can be i d e n t i f i e d i n the flow domain that are u s e f u l i n e f f i c i e n t design of sampling programs. A l t e r n a t i v e l y , can we e l i m i n a t e some p o t e n t i a l sampling l o c a t i o n s because t h e i r worth i s l e s s than the co s t of measurement? 7. I n v e s t i g a t e whether q u a l i t a t i v e g u i d e l i n e s f o r the design of piezometer networks f o r p r e d i c t i n g groundwater d i s c h a r g e can be developed. INTRODUCTION / 21 Use t h i s framework to compare d i f f e r e n t sampling programs and choose the best design among those c o n s i d e r e d . Test t h i s method on a r a p i d i n f i l t r a t i o n problem f o r a h y p o t h e t i c a l f i e l d s i t e . METHODOLOGY A number of d e c i s i o n s have to be made when d e s i g n i n g a sampling program. For the groundwater problem posed i n Fi g u r e 1, these d e c i s i o n s i n c l u d e the number of boreholes and t h e i r l o c a t i o n and the number of h y d r a u l i c c o n d u c t i v i t y measurements and t h e i r l o c a t i o n w i t h i n each borehole. There i s an i n f i n i t e set of sampling s t r a t e g i e s from which to choose. To a i d in the s e l e c t i o n of a p r e f e r r e d s t r a t e g y from a f i n i t e set of a l t e r n a t i v e s , an o b j e c t i v e f u n c t i o n can be used. In t h i s study, the o b j e c t i v e f u n c t i o n has two components; the cost of data c o l l e c t i o n and the worth of data. The cost of data c o l l e c t i o n , C^ ,, i s a l i n e a r f u n c t i o n of the d e c i s i o n v a r i a b l e s CS = f(^i>^i>> f ° r a given x[ where /V- i n t e g e r d e c i s i o n s such as the number of boreholes, the number of samples; Mi d e c i s i o n v a r i a b l e s such as meters d r i l l e d , meters cased; x, l o c a t i o n of the samples. 22 METHODOLOGY / 23 The Appendix c o n t a i n s a t a b l e of n o t a t i o n . For t h i s study, h y d r a u l i c c o n d u c t i v i t y i s measured at s e l e c t e d l o c a t i o n s on h y p o t h e t i c a l s i t e s f o r the purpose of p r e d i c t i n g volume flow through heterogeneous porous media. Data worth i s determined by c a l c u l a t i n g the economic l o s s , C^, a s s o c i a t e d with a poor p r e d i c t i o n of di s c h a r g e CL = f(Q,Qt) where Q estimated flow based on measurements and flow system a n a l y s i s ; Qt a c t u a l flow from h y p o t h e t i c a l f i e l d s i t e . The o b j e c t i v e of the sampling program i s to minimize the combined cost of sampling and the economic l o s s . An o b j e c t i v e f u n c t i o n e x p r e s s i n g t h i s concept can be w r i t t e n Z = f(Q,Q,,N.,Mi , xt) The estimated flow i s a f u n c t i o n of s i t e c h a r a c t e r i s t i c s and the number and l o c a t i o n of measurements For heterogeneous media, t h i s f u n c t i o n i s unknown. Because the r e l a t i o n s h i p expressed i n t h i s equation i s h i g h l y n o n - l i n e a r , a f i n i t e element model that preserves METHODOLOGY / 24 measurements i s used to compute estimates of flow. Values of Q are computed for d i f f e r e n t sampling programs, each based on a unique set of measurements (Afj , A4-, x-) . The sampling s t r a t e g y that y i e l d s the minimum value of Z i s the p r e f e r r e d a l t e r n a t i v e . By comparing d i f f e r e n t s t r a t e g i e s , g e n e r a l g u i d e l i n e s for the design of sampling programs can be in v e s t i g a t e d . T h i s method i s a p p l i e d to the problem of d e s i g n i n g a sampling program that measures h y d r a u l i c c o n d u c t i v i t y f o r p r e d i c t i n g discharge from a r a p i d i n f i l t r a t i o n pond. Rapid i n f i l t r a t i o n systems are u n l i n e d ponds that r e c e i v e m u n i c i p a l , i n d u s t r i a l or domestic wastewater [Swaney and Jackson, 1983; Ostendorf, 1986], As the e f f l u e n t moves through the subsurface, a v a r i e t y of chemical, p h y s i c a l , or b i o l o g i c a l processes can adsorb, break down, or p r e c i p i t a t e contaminants. The t r e a t e d e f f l u e n t may e v e n t u a l l y d i s c h a r g e i n t o s u r f a c e waters. In t h i s study, the h y p o t h e t i c a l v e r t i c a l c r o s s s e c t i o n shown i n Figu r e 2 i s used. The bottom boundary i s impermeable, r e p r e s e n t i n g the cont a c t between a more permeable u n i t above and a l e s s permeable u n i t , such as bedrock or an a q u i t a r d , below. The l e f t s i d e i s an impermeable boundary and rep r e s e n t s a g e o l o g i c c o n t a c t . The boundary on the r i g h t s i d e i s a r i v e r along which the h y d r a u l i c head i s co n s t a n t . METHODOLOGY / 25 Pond Ground surface 60 < > ID f Piezometer ports 32 64 96 128 160 192 224 256 288 320 HORIZONTAL DISTANCE, metres F i g . 2. P o s i t i o n of water t a b l e before and a f t e r pond c o n s t r u c t i o n . Suppose a d i s p o s a l pond of given dimension i s to be b u i l t at the ground s u r f a c e (dashed l i n e ) and we wish to design a sampling program to a i d i n p r e d i c t i n g the steady s t a t e c a p a c i t y of the pond. Because Qt i s known f o r a h y p o t h e t i c a l s i t e , our p r e d i c t i o n s can be compared to the a c t u a l value, Q(. The sampling phase w i l l i n v o l v e measuring h y d r a u l i c c o n d u c t i v i t y from a set of point samples, rather than o b t a i n i n g an average value from a pump t e s t . I t i s assumed that during t h i s pre-pond sampling phase, the i n i t i a l water t a b l e i s almost f l a t . Above the i n i t i a l water t a b l e , unsaturated samples are c o l l e c t e d , and below, s a t u r a t e d h y d r a u l i c c o n d u c t i v i t y measurements are made. In a l l cases, i t i s assumed that s u f f i c i e n t volumes of METHODOLOGY / 26 e f f l u e n t are s u p p l i e d to ensure that the pond i s f i l l e d with water at steady s t a t e and pr o v i d e s a constant head along the base of the pond [Wit her spoon and Narasimhan, 1972]. The constant head at the r i g h t hand sid e boundary i s lower, c r e a t i n g the head g r a d i e n t that d r i v e s flow. I t i s the s a t u r a t e d steady flow through the pond that we wish to p r e d i c t . The task i s to measure h y d r a u l i c c o n d u c t i v i t y at l o c a t i o n s that provide the g r e a t e s t r e d u c t i o n i n p r e d i c t i o n u n c e r t a i n t y of pond di s c h a r g e f o r the minimum c o s t . The importance of c h a r a c t e r i z i n g the s p a t i a l arrangement of the h e t e r o g e n e i t i e s i s ev a l u a t e d through c o n d i t i o n a l s i m u l a t i o n s . A h y p o t h e t i c a l , heterogeneous f i e l d of h y d r a u l i c c o n d u c t i v i t y v a l u e s i s generated and i s r e f e r r e d to as a referenc e case. A small set of p o t e n t i a l sampling s t r a t e g i e s i s s e l e c t e d . These sampling s t r a t e g i e s are o u t l i n e d l a t e r . The procedure f o r t e s t i n g the e f f e c t i v e n e s s of each s t r a t e g y i s as f o l l o w s . 1. Decide on the number and l o c a t i o n of boreholes and the number and l o c a t i o n of h y d r a u l i c c o n d u c t i v i t y measurements w i t h i n each borehole. 2. Use a s t o c h a s t i c model that p r e s e r v e s the measurements and i n c o r p o r a t e s the u n c e r t a i n t y where h y d r a u l i c c o n d u c t i v i t y i s unknown to generate a l a r g e number of d i f f e r e n t r e a l i z a t i o n s c o n s i s t e n t with the measurements. METHODOLOGY / 27 3. Solve each r e a l i z a t i o n d e t e r m i n i s t i c a l l y using a f i n i t e element model and a stream f u n c t i o n approach. 4. For each r e a l i z a t i o n c a l c u l a t e the pond d i s c h a r g e from the stream f u n c t i o n s o l u t i o n . 5. C a l c u l a t e the value of the l o s s f u n c t i o n by comparing the estimated d i s c h a r g e from the pond with the a c t u a l d i s c h a r g e from the r e f e r e n c e case. 6. C a l c u l a t e the cost of c o l l e c t i n g data. 7. C a l c u l a t e the value of the o b j e c t i v e f u n c t i o n , Z, by adding 5 and 6 above. The sampling s t r a t e g y which minimizes the o b j e c t i v e f u n c t i o n i s the p r e f e r r e d a l t e r n a t i v e . To t e s t the s e n s i t i v i t y of the s e l e c t e d a l t e r n a t i v e , the above procedure i s repeated f o r d i f f e r e n t r e f e r e n c e cases. STOCHASTIC MODEL S p a t i a l v a r i a b i l i t y A s t o c h a s t i c model i s used to i n c o r p o r a t e h e t e r o g e n e i t y and u n c e r t a i n t y i n the d i s t r i b u t i o n of h y d r a u l i c c o n d u c t i v i t y . The model assumes the porous medium i s s t a t i s t i c a l l y homogeneous and uses three parameters to c h a r a c t e r i z e s p a t i a l v a r i a b i l i t y . These are the mean, va r i a n c e and c o v a r i a n c e . The assumption of s t a t i s t i c a l homogeneity im p l i e s that the expected value of h y d r a u l i c METHODOLOGY / 28 c o n d u c t i v i t y and the v a r i a n c e are constant i n space. In a d d i t i o n , the co v a r i a n c e depends only on the d i s t a n c e s e p a r a t i n g two p o i n t s and not on t h e i r o r i e n t a t i o n or p o s i t i o n in space. The s t o c h a s t i c process from which each h y d r a u l i c c o n d u c t i v i t y r e a l i z a t i o n i s d e r i v e d , i s assumed to be e r g o d i c . T h i s i m p l i e s that the ensemble s t a t i s t i c s can be determined by d e t a i l e d sampling w i t h i n any one r e a l i z a t i o n , s p e c i f i c a l l y , the one r e a l i z a t i o n that i s a v a i l a b l e to us i n a r e a l case. Many r e s e a r c h e r s have observed that h y d r a u l i c c o n d u c t i v i t y data can be d e s c r i b e d by a lognormal d i s t r i b u t i o n [ F r e e z e , 1975; Smith, 1981]. In order to work with a normally d i s t r i b u t e d parameter, we use the transform Y - l o g 1 0 K. The v a r i a n c e i n Y i s a measure of the magnitude of v a r i a b i l i t y i n h y d r a u l i c c o n d u c t i v i t y . I t has been found to range between 0.2, f o r r e l a t i v e l y homogeneous g e o l o g i c u n i t s , to 1.6, f o r heterogeneous media [ F r e e z e , 1975; El-Kadi and Brutsaert, 1985]. The s p a t i a l c o n t i n u i t y i s expressed using an autocova r i a n c e f u n c t i o n . In t h i s study the autocovariance f u n c t i o n i s assumed to be an e x p o n e n t i a l f u n c t i o n c o v ( l x , l z ) = a2exp{-V[U x/\x)2 + (/ zA z) 2]) ( 1 ) where /;- i s the l a g or s e p a r a t i o n d i s t a n c e i n d i r e c t i o n / between two p o i n t s . The i n t e g r a l s c a l e , X,, i s a measure of METHODOLOGY / 29 the average d i s t a n c e over which h y d r a u l i c c o n d u c t i v i t y v a l u e s are c o r r e l a t e d i n d i r e c t i o n / [Bakr et al., 1978]. I t i s d e f i n e d by the d i s t a n c e at which the c o r r e l a t i o n (cov/o2) decays to e' 1 or 0.3679 [Gel har, 1986], The l a r g e r the i n t e g r a l s c a l e , the gr e a t e r the d i s t a n c e over which s p a t i a l c o n t i n u i t y i n h y d r a u l i c c o n d u c t i v i t y i s observed. The a n i s o t r o p y inherent i n l a y e r e d d e p o s i t s i s expressed by s p e c i f y i n g an i n t e g r a l s c a l e that i s l a r g e r p a r a l l e l to l a y e r s than p e r p e n d i c u l a r to l a y e r s . The i n t e g r a l s c a l e i s small compared to the s i z e of the flow domain. If h e t e r o g e n e i t i e s e x i s t on a s c a l e s i m i l a r to the s c a l e of the flow domain, they must be modeled as d i s t i n c t u n i t s with separate s t a t i s t i c a l p r o p e r t i e s [Smith, 1984]. G e n e r a t i n g the r e f e r e n c e case There are s e v e r a l ways of gene r a t i n g s p a t i a l l y c o r r e l a t e d random f i e l d s (nearest neighbor, matrix method, t u r n i n g bands). In t h i s study, the matrix method of Scheuer and Stoller [1962] and Clifton and Neuman [1982] i s used to generate h y d r a u l i c c o n d u c t i v i t y f i e l d s . The region being modeled i s d i v i d e d i n t o r e c t a n g u l a r b l o c k s of uniform s i z e (Figure 3 ) . The r e c t a n g u l a r blocks are longer i n the h o r i z o n t a l d i r e c t i o n , r e f l e c t i n g the inherent a n i s o t r o p y of l a y e r e d d e p o s i t s . A pxp c o v a r i a n c e matrix, V, i s set up where p i s the number of s t o c h a s t i c METHODOLOGY / 30 60 30 0 >* > L 1 T H 0 32 64 96 128 160 192 224 256 288 320 HORIZONTAL DISTANCE, metres F i g . 3. S t o c h a s t i c b l o c k s . b l o c k s . The matrix i s formed using the autocovariance f u n c t i o n ( 1 ) . The entry V-t • r e f e r s to the cov a r i a n c e between h y d r a u l i c c o n d u c t i v i t y values i n block / and block j . The matrix i s symmetric with o2 on the d i a g o n a l . We want to f i n d a matrix M such that MMT = V. A Cholesky decomposition of V y i e l d s M. If we m u l t i p l y M by a random vect o r R from the N(0 , 1 ) d i s t r i b u t i o n we have.a random f i e l d with a N(0,o 2) d i s t r i b u t i o n . To convert to a N(y,a 2) d i s t r i b u t i o n we have only to add the d e s i r e d mean, u. These r e a l i z a t i o n s form the h y p o t h e t i c a l r e f e r e n c e cases and are u n c o n d i t i o n a l because valu e s are not f i x e d i n any l o c a t i o n . The average value of log K i n a l l blocks forming the referenc e case, Y , and standard d e v i a t i o n , Sr, d i f f e r s l i g h t l y from the ensemble s t a t i s t i c s , u and a, from which they were d e r i v e d due to s t a t i s t i c a l f l u c t u a t i o n s . METHODOLOGY / 31 C o n d i t i o n a l s i m u l a t i o n H y d r a u l i c c o n d u c t i v i t y v a l u e s from the re f e r e n c e case are sampled and the sample average, Y' , and standard d e v i a t i o n , S™, c a l c u l a t e d . Except by chance, these s t a t i s t i c s w i l l be d i f f e r e n t from the re f e r e n c e case Yr and Sry due to s t a t i s t i c a l v a r i a t i o n s . The measurements of h y d r a u l i c c o n d u c t i v i t y are assumed to be f r e e from e r r o r . C o n d i t i o n a l r e a l i z a t i o n s are generated that preserve both the measurements and t h e i r average value, Y , and a s s i g n random valu e s at unknown b l o c k s . When the number of measurements i s s m a l l , a ccurate estimates of the standard d e v i a t i o n of the re f e r e n c e case are not p o s s i b l e . For t h i s reason a d e s i r e d standard d e v i a t i o n f o r the r e a l i z a t i o n can be s p e c i f i e d , such as Sr i n s t e a d of S™. C o n d i t i o n a l r e a l i z a t i o n s use k r i g i n g , which i s a type of i n t e r p o l a t i o n , to extend the i n f l u e n c e of measured val u e s to the surrounding areas [Journel and Huijbregis, 1978; Delhomme, 1979; Clifton and Neuman, 1982]. The u n c e r t a i n t y i n h y d r a u l i c c o n d u c t i v i t y decreases the c l o s e r a block i s to a measurement. Because of t h i s , c o n d i t i o n a l models are no n s t a t i o n a r y . The f o l l o w i n g d e s c r i p t i o n of the technique fo r forming a c o n d i t i o n a l r e a l i z a t i o n i s based on the work of Smith and Schwartz [19816]. METHODOLOGY / 32 The k r i g i n g estimate i s d e f i n e d as Y*(i,j) = l\tYt + Ym (2) where / row number of s t o c h a s t i c block; j column number of s t o c h a s t i c block; n number of measurements; X. k r i g i n g weights. Ym mean of measured Y val u e s ; Yi measured l o g K v a l u e s a d j u s t e d t o a mean of zero; The k r i g i n g weights f o r each block, X> , are found by s o l v i n g the system of equations Au = E A = ClYuY,) C(YUY2) C(y a ,F i ) C(Y2,Y2) C{YniYx) C(Yn,Y2) 1 1 C(YuYn) V C(Y2,Yn) 1 C(Yn,Yn) 1 1 0 u = / A , \ An V A* J and {C(Y1,Y*(i,j)\ C(Y2,Y*(i,j) E = C(Yn,Y'(i,j) 1 METHODOLOGY / 33 where u Lagrange m u l t i p l i e r ; C(YX,Y2) c o v a r i a n c e between data p o i n t 1 and data point 2; C(Y,,Y*(i, j)) c o v a r i a n c e between data p o i n t 1 and the k r i g i n g estimate f o r unknown block An u n c o n d i t i o n a l r e a l i z a t i o n , s(i,j), i s formed using e i t h e r the sample v a r i a n c e , (S™)2, or the po p u l a t i o n v a r i a n c e , oy2. The k r i g i n g estimate from the u n c o n d i t i o n a l r e a l i z a t i o n , s i s formed by k r i g i n g the u n c o n d i t i o n a l values at the n sampling l o c a t i o n s using the k r i g i n g weights found i n ( 2 ) . A c o n d i t i o n a l r e a l i z a t i o n , Y (i,j), i s then formed by generating sUU ,j ), s*(i,j) and s o l v i n g f o r Y c ( i , j ) = Y*(i,j) + [ s j i . j ) - s*(i,j)] (3) where it Y (i,j) k r i g i n g estimate of l o g K f o r block i,j formed by i n t e r p o l a t i n g between measured values; su(i>j) u n c o n d i t i o n a l r e a l i z a t i o n of l o g K; METHODOLOGY / 34 s (i,j) k r i g e d " u n c o n d i t i o n a l " r e a l i z a t i o n formed by k r i g i n g the l o g K values from the u n c o n d i t i o n a l s i m u l a t i o n at the measurement l o c a t i o n s . In the Monte C a r l o technique, a c o n d i t i o n a l s i m u l a t i o n i s formed from a l a r g e number of c o n d i t i o n a l r e a l i z a t i o n s i n which h y d r a u l i c c o n d u c t i v i t y values remain constant f o r measured blocks but vary f o r a l l other b l o c k s . As i t i s expensive to decompose l a r g e m a t r i c e s , the matrix, M, i s found only f o r the f i r s t u n c o n d i t i o n a l r e a l i z a t i o n and the r e s u l t s t o r e d . One only needs to generate a d i f f e r e n t random ve c t o r R and m u l t i p l y by M f o r each subsequent u n c o n d i t i o n a l r e a l i z a t i o n . The normally d i s t r i b u t e d Y values from each c o n d i t i o n a l r e a l i z a t i o n are transformed i n t o lognormally d i s t r i b u t e d K values using the e x p o n e n t i a l transform K(i,j) = exp[2.3026 Y(i,j)]. These values of h y d r a u l i c c o n d u c t i v i t y are used i n the f i n i t e element model to solve f o r steady s t a t e flow. METHODOLOGY / 35 SOLUTION Stream f u n c t i o n s A c o n d i t i o n a l r e a l i z a t i o n generates a f i e l d of h y d r a u l i c c o n d u c t i v i t y values that can be used d e t e r m i n i s t i c a l l y i n a f i n i t e element model that s o l v e s f o r stream f u n c t i o n s . Stream f u n c t i o n s lend themselves q u i t e r e a d i l y to determining discharge a c r o s s boundaries of the flow domain [Frind and Mat a n g a, 1 98 5 j Frind el a I, 1985]. For each r e a l i z a t i o n , the estimate of pond d i s c h a r g e i s st o r e d along with the average value and standard d e v i a t i o n of the c o n d i t i o n a l h y d r a u l i c c o n d u c t i v i t y f i e l d . These s t a t i s t i c s are accumulated over the l a r g e number of r e a l i z a t i o n s that form a s i m u l a t i o n and then the average estimate of pond d i s c h a r g e and standard d e v i a t i o n f o r the s i m u l a t i o n are c a l c u l a t e d . There are a number of advantages to using stream f u n c t i o n s as an a l t e r n a t i v e to s o l v i n g f o r h y d r a u l i c heads. It i s e a s i e r to v i s u a l i z e p a t t e r n s of flow from maps of stream e q u i p o t e n t i a l s because these c o i n c i d e with s t r e a m l i n e s . Boundary f l u x e s are more e a s i l y and more a c c u r a t e l y c a l c u l a t e d d i r e c t l y from the s o l u t i o n because values of stream f u n c t i o n s at boundary nodes represent cumulative f l u x e s . Frind and Mat anga [1985] d e s c r i b e the f i n i t e element f o r m u l a t i o n f o r stream f u n c t i o n s i n d e t a i l . A METHODOLOGY / 36 b r i e f summary of the d e r i v a t i o n of the governing equations f o l l o w s . A stream f u n c t i o n i s d e f i n e d as ^ = \js(x,z) ' with dimensions of volume per u n i t time per un i t l e n g t h [L2/T]. Assuming the co o r d i n a t e axes are p a r a l l e l to the p r i n c i p a l d i r e c t i o n s of h y d r a u l i c c o n d u c t i v i t y , the equation l i n k i n g stream p o t e n t i a l s to h y d r a u l i c p o t e n t i a l i s JC d\ / dtp dz J \ dx where qi i s f l u x [L/T] i n d i r e c t i o n /, K(. i s an element of the h y d r a u l i c c o n d u c t i v i t y tensor [L/T] and

= K- 1 q (7) METHODOLOGY / 37 S u b s t i t u t i n g (7) and (5) i n t o (6), expanding, and then s u b s t i t u t i n g (4), we have the governing equation 3 1 3 ^ 9 1 — [ ] + — [ ] = 0 (8) dx AT bx dz Kvv 3z Z Z XX where 1/^// i s h y d r a u l i c r e s i s t i v i t y . Boundary c o n d i t i o n s Boundaries along which the value of the stream f u n c t i o n i s known are Type I (or D i r i c h l e t ) boundaries. Type II (or Neumann) boundary c o n d i t i o n s occur where the normal g r a d i e n t in the stream f u n c t i o n i s s p e c i f i e d . T h i s second type boundary i s expressed by the equation n«VuV = -r«V0 where n and r are the normal and t a n g e n t i a l u n i t v e c t o r s . The normal component of the g r a d i e n t i n the stream f u n c t i o n can be found by the rate of change i n heads tangent to the boundary. Along constant head boundaries, n«Vu/ = 0. Di scharge Stream f u n c t i o n s can be v i s u a l i z e d as a map of volume dis c h a r g e through the flow domain. The flow or Darcy d i s c h a r g e , A£, through a stream tube i s c a l c u l a t e d as the d i f f e r e n c e between the value of the stream f u n c t i o n s on the bounding s t r e a m l i n e s . The boundary f l u x or s p e c i f i c d i s c h arge [L/T] along a boundary segment, Ar, of a stream METHODOLOGY / 38 tube, AuV, i s qb = W/L\T (9) The s o l u t i o n technique f o l l o w s the G a l e r k i n procedure with l i n e a r b a s i s f u n c t i o n s [Frind and Matanga, 1985]. Equation (9) i s used to sol v e f o r flow through the pond. E V A L U A T I O N OF SAMPLING S T R A T E G I E S Data c o l l e c t i o n costs Once the number and l o c a t i o n of measurements f o r each sampling s t r a t e g y i s s e l e c t e d , the cost of c o l l e c t i n g data can be determined. The co s t i s a f u n c t i o n of the d e c i s i o n v a r i a b l e s and Mi and co s t c o e f f i c i e n t s Cz- . The cost of c o l l e c t i n g h y d r a u l i c c o n d u c t i v i t y data i s C M U A + c r u , + c „ N r , + C J „ + C.N, + C, a a c c p p u u n n I m (10) where d r i l l i n g c o s t s ($/m); cost of piezometer casing($/m); cost f o r each s a t u r a t e d sample (screens,sampling, zone development); C u cost f o r each unsaturated zone sample ( c o l l e c t i n g , s h i p p i n g , t e s t i n g ) ; METHODOLOGY / 39 C^ c o s t per borehole ( c o l l a r completions, expenses and t r a v e l f o r f i e l d h y d r o g e o l o g i s t ) ; Cm m i s c e l l a n e o u s c o s t s ( i . e . d r i l l i n g m o b i l i z a t i o n / d e m o b i l i z a t i o n , compressor r e n t a l ) ; Mj meters d r i l l e d ; Mc meters of c a s i n g ; N number of piezometers; N number of unsaturated samples; N^ number of boreholes; R e c a l l that d u r i n g the pre-pond sampling phase, there i s a la r g e unsaturated zone. Core samples would be c o l l e c t e d from the unsaturated zone i n the pre-development stage f o r determining h y d r a u l i c c o n d u c t i v i t y u s i n g , g r a i n s i z e a n a l y s i s or permeameter t e s t s . In the s a t u r a t e d zone, s l u g or b a i l t e s t s would be used f o r measuring h y d r a u l i c c o n d u c t i v i t y . The c o s t s assume that m u l t i p l e piezometers would be bundled together i n a s i n g l e borehole. METHODOLOGY / 40 Data worth A l o s s f u n c t i o n i s used to determine the economic worth of data that are c o l l e c t e d . T h i s f u n c t i o n expresses the economic l o s s i n c u r r e d from poor estimates of pond d i s c h a r g e . I f a sampling s t r a t e g y r e s u l t s i n a good estimate of d i s c h a r g e f o r a given design, the pond w i l l perform as planned and the l o s s f u n c t i o n w i l l be minimal. The l o s s f u n c t i o n would only be zero i f the estimate of di s c h a r g e f o r each r e a l i z a t i o n i n a c o n d i t i o n a l s i m u l a t i o n i s equal to the a c t u a l d i s c h a r g e . I f , however, the sampling s t r a t e g y r e s u l t s i n poor estimates of dis c h a r g e f o r a given design, the pond w i l l perform p o o r l y and r e s u l t in economic l o s s e s . The more the estimates of pond d i s c h a r g e d i f f e r from the tr u e value, the g r e a t e r the economic l o s s . Using t h i s framework, the worth of data can be c a l c u l a t e d and compared to the c o s t of c o l l e c t i n g these data. The root mean square e r r o r has been chosen i n t h i s study as the c r i t e r i o n by which to judge the worth of data. In the f o l l o w i n g d i s c u s s i o n , the root mean square e r r o r i s d e r i v e d by f i n d i n g the square root of the expected value of the e s t i m a t i o n e r r o r squared. Let us begin by c o n s i d e r i n g the square of the e s t i m a t i o n e r r o r f o r s i n g l e r e a l i z a t i o n s , which i s the d e v i a t i o n of estimates of pond d i s c h a r g e , Q, from the true v a l u e , Qt (Figure 4 ) . T h i s p a r a b o l i c f u n c t i o n i s expressed METHODOLOGY / 41 Underestimate Overestimate F i g . 4. Function e x p r e s s i n g the square of the b i a s . as f(Q) = (Q ~ Q,)2 (11) where Q estimate of pond d i s c h a r g e from a c o n d i t i o n a l r e a l i z a t i o n Qt a c t u a l d i s c h a r g e from the pond f o r a ref e r e n c e case. METHODOLOGY / 42 For the set of r e a l i z a t i o n s forming a s i m u l a t i o n , we have E[f(Q) ] = ;(Q - Qt )2 P(Q) dQ = 5Q2P(Q)dQ-2Qt jQP(Q) dQ+Qt 2jP(Q)dQ R e c a l l i n g $P(Q)dQ = / Q = SQP(Q)dQ S q 2 = / (G " Q)2P(Q)dQ = lQ2P(Q)dQ - Q2 we have E[f(Q)] = Q2 + S 2 - 2QQ{ + Qt 2 E[f(Q)] = (Q - Q. )2 + S 2 i q where Q average value of pond d i s c h a r g e f o r a s imulat i o n . T h i s equation i s the mean square e r r o r [Benjamin and Cornell, 1970]. The root mean square e r r o r (RMSE) i s RMSE = •[(Q - Qt )2 + S g 2 ] (12) The l o s s f u n c t i o n used i n t h i s study i s CL = C[V/[ (Q ~ Q t ) 2 + S 2] (13) METHODOLOGY / 43 where cost c o e f f i c i e n t based on the c o s t of a l t e r n a t i v e treatment, the s i z e of the pond, and the time h o r i z o n the pond i s designed f o r . T h i s equation expresses the economic l o s s i n terms of the root mean square e r r o r . The b i a s , or d e v i a t i o n of the average estimate from the tr u e v a l u e , and the v a r i a b i l i t y of the estimates are accounted f o r . The cost i n c r e a s e s as the square root of the sum of the square of the d e v i a t i o n , i n d i c a t i n g that the f u r t h e r the estimates are from the tr u e v a l u e , the gr e a t e r the economic l o s s . In the event that the average estimate of d i s c h a r g e i s equal to Qt , the expected c o s t w i l l r e f l e c t the v a r i a b i l i t y i n i n d i v i d u a l estimates of di s c h a r g e f o r each r e a l i z a t i o n . The c o e f f i c i e n t of v a r i a t i o n , v = a/u, i s f r e q u e n t l y used as an i n d i c a t i o n of p r e d i c t i o n u n c e r t a i n t y [Smith and Schwartz, 1981a; Smith and Schwartz, 19816]. For t h i s study, RMSE appears to be a b e t t e r c r i t e r i o n f o r determining the accuracy of di s c h a r g e p r e d i c t i o n s than v = S /Q. The root mean square e r r o r weights both the b i a s and 5 while the c o e f f i c i e n t of v a r i a t i o n ignores the b i a s and i s i n f l u e n c e d by S and Q. The c o e f f i c i e n t of v a r i a t i o n tends to favor o v e r estimates of Qt . I f the average estimate of d i s c h a r g e , Q, i s an overestimate, the c o e f f i c i e n t of v a r i a t i o n w i l l be METHODOLOGY / 44 • lower than v f o r a s i m u l a t i o n that underestimates Qf but q « has the same magnitude of b i a s and same S^. The RMSE, however, would be the same f o r both s i m u l a t i o n s . Another problem with the c o e f f i c i e n t of v a r i a t i o n i s that a sampling scheme that overestimates Q. would have a lower v than a _ q sampling scheme that has no b i a s (Q = Q.), but the same S . i q For these reasons, RMSE i s used in t h i s study as the c r i t e r i o n of p r e d i c t i o n u n c e r t a i n t y . It i s not the purpose of t h i s study to determine an exact form f o r the l o s s f u n c t i o n , but ra t h e r to i l l u s t r a t e the concept. The exact form of the f u n c t i o n would be determined on a case by case b a s i s and depends on such f a c t o r s as s i t e l o c a t i o n , type of l i q u i d d i s posed of and cost of a l t e r n a t i v e treatments. The l o s s f u n c t i o n need not be symmetric. A symmetric f u n c t i o n has been chosen here to ensure that the r e s u l t s are not b i a s e d toward e i t h e r o v e r e s t i m a t i n g or underestimating the pond c a p a c i t y . In t h i s study, we assume the dimensions of the d i s p o s a l pond are f i x e d . The purpose of the sampling program i s to p r e d i c t d i s c h a r g e from the pond, and to use t h i s i n f o r m a t i o n i n the design of a d i s p o s a l f a c i l i t y . From the sampling program, flow through the pond i s p r e d i c t e d . On the b a s i s of p r e d i c t e d d i s c h a r g e , d e c i s i o n s would be made on the s i z e of f a c i l i t y to b u i l d or the volume of waste to accept or produce at the f a c i l i t y . These d e c i s i o n s c o u l d r e s u l t i n METHODOLOGY / 45 economic l o s s i f the p r e d i c t i o n has l a r g e e r r o r . In the case of an overestimate of pond c a p a c i t y , l e s s water i n f i l t r a t e s through the base of the pond than p r e d i c t e d . If i t i s necessary to dispose of a g r e a t e r volume of e f f l u e n t than the pond i s p r e d i c t e d to handle, then a l t e r n a t i v e d i s p o s a l methods must be found or the pond has to be expanded. E i t h e r of these a l t e r n a t i v e s r e s u l t i n a d d i t i o n a l expenses, and are manifested as an i n c r e a s e i n the l o s s f u n c t i o n . Suppose that s i t e e x p l o r a t i o n leads to underestimates of pond d i s c h a r g e . In t h i s case, more water i n f i l t r a t e s through the base of the pond than p r e d i c t e d . S e v e r a l s c e n a r i o s c o u l d r e s u l t i n economic l o s s . I f the s i z e of the f a c i l i t y i s l i m i t e d and the maximum throughput exceeds the c a p a c i t y of the f a c i l i t y , then more money would have been spent i n c o n s t r u c t i n g the pond than necessary. Or perhaps the owner-operator has c o n t r a c t e d to accept a f i x e d volume of waste when, in f a c t , a l a r g e r volume c o u l d be handled. Another s c e n a r i o that r e s u l t s i n economic l o s s e s i s i f the ra t e of i n f i l t r a t i o n from the pond i s high and the e f f l u e n t moves so r a p i d l y through the subsurface that i t i s improperly t r e a t e d . I t would then have to be c o l l e c t e d and r e t r e a t e d . Another p o s s i b i l i t y f o r economic l o s s i s i f remedial a c t i o n i s r e q u i r e d because of slope f a i l u r e i n di s c h a r g e areas due to higher than p r e d i c t e d v e l o c i t i e s and g r a d i e n t s . T h i s may a l s o r e q u i r e lowering the pond l e v e l METHODOLOGY / 46 while the water t a b l e mound d i s s i p a t e s . A l l of these examples would r e s u l t i n economic l o s s e s from poor e x p l o r a t i o n s t r a t e g i e s . O b j e c t i v e f u n c t i o n The c o s t of the sampling program i s combined with the l o s s f u n c t i o n in an o b j e c t i v e f u n c t i o n mi n Z = Cs + C / V/[ (Q - Qt )2 + S q 2 ] (14) M i n i m i z a t i o n of t h i s o b j e c t i v e f u n c t i o n p r o v i d e s a simple method of comparing d i f f e r e n t sampling schemes to determine a p r e f e r r e d a l t e r n a t i v e . In a d d i t i o n t h i s framework compares the worth of data to the cost of data c o l l e c t i o n . SOURCE OF COMPUTER CODES The computer model used i n t h i s study, i s a combination of r o u t i n e s developed by the author and e x i s t i n g computer codes. The s t o c h a s t i c flow model developed by Smith and Schwartz [19816] was m o d i f i e d s u b s t a n t i a l l y . The r o u t i n e f o r ge n e r a t i n g s p a t i a l l y c o r r e l a t e d random f i e l d s was a c q u i r e d by J o e l Massmann while at the U n i v e r s i t y of A r i z o n a and made a v a i l a b l e to the author. The stream f u n c t i o n s o l u t i o n was adapted from a subroutine w r i t t e n by L e s l i e Smith. A f r e e s u r f a c e model, developed by C r a i g F o r s t e r , was used to f i n d the water t a b l e f o r a homogeneous case and to generate the METHODOLOGY / 47 mesh. The a l g o r i t h m f o r c a l c u l a t i n g the co s t of c o l l e c t i n g data and the value of the l o s s f u n c t i o n was developed by the author. A l l the a c q u i r e d r o u t i n e s were m o d i f i e d , adapted to the FPS-164 atta c h e d processor at the U n i v e r s i t y of B r i t i s h Columbia, and v a l i d a t e d by the author. COMPUTER SIMULATIONS PHYSICAL MODEL Boundary value problem F i g u r e 5 shows the flow domain used f o r a l l s i m u l a t i o n s . The upstream p o r t i o n of the upper boundary r e p r e s e n t s an i n f i l t r a t i o n or d i s p o s a l pond of f i x e d l e n g t h . The remainder of the upper boundary i s a water t a b l e . The p o s i t i o n of the water t a b l e i s h e l d constant through a l l the s i m u l a t i o n s . I t was found using a f r e e s u r f a c e model f o r homogeneous cases with a h y d r a u l i c c o n d u c t i v i t y value between 10" 1 m/s and 10" 5 m/s. Because of the unique arrangement of the h y d r a u l i c c o n d u c t i v i t y f i e l d i n each r e a l i z a t i o n , one expects the water t a b l e p o s i t i o n to vary between r e a l i z a t i o n s . A s t o c h a s t i c f r e e s u r f a c e model, however, was beyond the scope of t h i s study. The water t a b l e c o n f i g u r a t i o n was f i x e d f o r a l l the r e a l i z a t i o n s forming a Monte C a r l o s i m u l a t i o n . The l e f t s i d e and bottom are impermeable boundaries along which the stream f u n c t i o n i s zero . Along the constant head boundary on the r i g h t s i d e , the normal g r a d i e n t i n the stream f u n c t i o n i s zero . The l e n g t h of the flow domain i s 320 meters with a 30 meter head drop between the pond and the r i g h t s i d e boundary. Equation 8 f o r steady s a t u r a t e d 48 COMPUTER SIMULATIONS / 49 Z L = 60m -y Z » 30m n-Vf - o Type II 320m F i g . 5. Groundwater flow domain and boundary con d i t i o n s . flow i n two dimensions i s used. F i n i t e element g r i d A f i n i t e element model i s used to solve each c o n d i t i o n a l r e a l i z a t i o n of h y d r a u l i c c o n d u c t i v i t y f o r the nodal v a l u e s of the stream f u n c t i o n . The f i n i t e element g r i d ( F i g u r e 6) i s denser where v e l o c i t i e s are expected to be h i g h e r . To ensure that each element l i e s e n t i r e l y w i t h i n a s i n g l e s t o c h a s t i c block, and to f u l f i l l i n t e r n a l requirements of the mesh generator, a d d i t i o n a l v e r t i c a l l i n e s are i n s e r t e d . In t h i s way, a value of h y d r a u l i c c o n d u c t i v i t y can be a s s i g n e d to each element. COMPUTER SIMULATIONS / 50 0 30 60 90 120 150 180 210 240 270 300 HORIZONTAL DISTANCE, metres F i g . 6. F i n i t e element g r i d . S T O C H A S T I C S I M U L A T I O N S Reference cases Table 1 shows the input values of u a X , and X y y z used to generate each r e f e r e n c e case. The f i r s t f i v e r e f e r e n c e cases (numbers 1 through 5) were generated from the same p r o b a b i l i t y model f o r h y d r a u l i c c o n d u c t i v i t y . The mean, a y l which i s -3.0 f o r K measured i n m/s {K = 10" 3 m/s), i s t y p i c a l of c l e a n sands [Freeze and Cherry, 1979]. R e c a l l i n g that measurements of oy commonly range from 0.2 to 1.6, a standard d e v i a t i o n of 0.2 was used to represent f a i r l y uniform media. To r e f l e c t d e p o s i t i o n a l l a y e r i n g , a h o r i z o n t a l i n t e g r a l s c a l e , X x, of 32 m and a v e r t i c a l i n t e g r a l s c a l e of 12m were chosen. These i n t e g r a l s c a l e s COMPUTER SIMULATIONS / 51 TABLE 1. Input and Output f o r Reference Cases Ref. Model Input M odel Output Case No. A, A* Yr ^y Qt 1 -3.0 0.2 32 12 -2.93 0.20 0.73 2 -3.0 0.2 32 12 -2.96 0.17 0.46 3 -3.0 0.2 32 12 -3.01 0.20 0.48 4 -3.0 0.2 32 12 -3.05 0.18 0.33 5 -3.0 0.2 32 12 -2.95 0.19 0.51 6 -3.0 0.5 32 12 -3.03 0.41 0.3S 7 -3.0 0.5 32 12 -2.67 0.44 1.45 8 -3.0 0.5 32 12 -2.90 0.47 0.55 9 -3.0 0.5 32 12 -3.04 0.46 0.41 10 -3.0 0.5 32 12 -3.08 0.41 0.42 11 -3.0 0.5 16 6 -2.92 0.55 0.93 12 -3.0 0.5 16 6 -3.09 0.46 0.82 13 -3.0 0.5 16 6 -3.11 0.48 0.47 14 -3.0 0.5 160 6 -2.89 0.51 0.62 15 -3.0 0.5 160 6 -3.15 0.43 0.37 16 -3.0 0.5 160 6 -3.06 0.49 0.67 17 -3.0 0.5 160 6 -2.06 0.47 0.71 * avi Yr, Sy for K in m/«; A in m; Qt in (m3/s x 10 2) represent strong s p a t i a l a u t o c o r r e l a t i o n f o r the s c a l e of the problem and c r e a t e reasonably l a r g e zones with s i m i l a r values of h y d r a u l i c c o n d u c t i v i t y . Reference cases 6 through 10 were designed to i n v e s t i g a t e the e f f e c t s of l a r g e r standard d e v i a t i o n . Input parameters were chosen that t y p i f y a sand with lenses of s i l t y sand and g r a v e l . Because of the l a r g e r standard d e v i a t i o n f o r these r e f e r e n c e cases (a^,=0.5) than f o r the prev i o u s r e f e r e n c e cases, the magnitude of the c o n t r a s t between h y d r a u l i c c o n d u c t i v i t y values i s g r e a t e r . COMPUTER SIMULATIONS / 52 Reference cases 11 through 13 have s m a l l e r i n t e g r a l s c a l e s (X =16m, X =6m) than the other r e f e r e n c e cases. They X z J t y p i f y a sandy a q u i f e r with s m a l l e r lenses of s i l t y sand and g r a v e l . In order to i n v e s t i g a t e the e f f e c t i v e n e s s of sampling s t r a t e g i e s i n s t a t i f i e d media, r e f e r e n c e cases 14 through 17 were generated using a l a r g e h o r i z o n t a l i n t e g r a l s c a l e (Xx=160m) and a small v e r t i c a l i n t e g r a l s c a l e (X z=6m). S p a t i a l averaging When a sampling s t r a t e g y i s a p p l i e d to a r e f e r e n c e case, a value f o r h y d r a u l i c c o n d u c t i v i t y i s obtained at each measurement p o i n t . In a s s i g n i n g p o i n t measurements of h y d r a u l i c c o n d u c t i v i t y from a f i e l d s i t e to b l o c k s , s p a t i a l averaging must be performed to preserve the s t a t i s t i c s of the p o i n t v a l u e s . In t h i s study we assume that h y d r a u l i c c o n d u c t i v i t y measurements apply to blocks d i r e c t l y . T h i s i s j u s t i f i e d because the flow f i e l d s of the r e f e r e n c e cases we are attempting to r e p l i c a t e are based on d i s c r e t e values of h y d r a u l i c c o n d u c t i v i t y r a t h e r than a continuum. Sample g r i d s The r e f e r e n c e cases are sampled using v a r i o u s s t r a t e g i e s . The purpose of the s t r a t e g i e s i s to i n v e s t i g a t e g u i d e l i n e s f o r measuring h y d r a u l i c c o n d u c t i v i t y f o r use i n p r e d i c t i n g d i s c h a r g e . S e v e r a l q u e s t i o n s that focus on how to COMPUTER SIMULATIONS / 53 measure h y d r a u l i c c o n d u c t i v i t y are addressed. 1. Is there a p r e f e r r e d l o c a t i o n f o r i n i t i a l boreholes? 2. Are there more important regions to sample? 3. Where should a second borehole be l o c a t e d given the l o c a t i o n of an i n i t i a l borehole? 4. In s t r a t i f i e d media, what are the t r a d e o f f s between c o l l e c t i n g more measurements i n boreholes and d r i l l i n g new boreholes? 5. How e f f e c t i v e are m u l t i p l e boreholes i n reducing p r e d i c t i o n u n c e r t a i n t y ? At what p o i n t does the in c r e a s e i n c o s t s of a d d i t i o n a l boreholes outweigh the i n f o r m a t i o n gain? In a d d i t i o n , the s e n s i t i v i t y of the r e s u l t s to changes i n input parameters i s i n v e s t i g a t e d . The f i r s t s e r i e s of s i m u l a t i o n s i s designed to determine whether a p r e f e r r e d l o c a t i o n e x i s t s f o r i n i t i a l samples. Before d r i l l i n g commences i n the f i e l d , the s i t e i n v e s t i g a t o r assumes there i s a s i n g l e g e o l o g i c u n i t . With l i t t l e p r e l i m i n a r y i n f o r m a t i o n a v a i l a b l e , a d e c i s i o n must be made where to begin d r i l l i n g . F i g u r e 7 shows four p o s s i b l e l o c a t i o n s f o r an i n i t i a l borehole. Twelve measurements of h y d r a u l i c c o n d u c t i v i t y , spaced three meters apa r t , are made down the borehole. Because the s t o c h a s t i c blocks are three COMPUTER SIMULATIONS / 54 HORIZONTAL DISTANCE, metres F i g . 7. L o c a t i o n of measurements i n i n i t i a l boreholes A through D. meters high, t h i s corresponds to one measurement per bl o c k . Borehole A i s l o c a t e d below the middle of the pond, B i s below the edge of the pond, C i s i n the middle of the region of flow and D i s near the constant head boundary. The second set of s i m u l a t i o n s i s designed to address q u e s t i o n 2 above. I n t u i t i v e l y , one may f e e l that measurements of h y d r a u l i c c o n d u c t i v i t y in one area may have a gr e a t e r e f f e c t on reducing u n c e r t a i n t y i n discharge estimates than measurements i n other areas. For example, we may b e l i e v e that Region F in F i g u r e 8, comprising the d e p o s i t s d i r e c t l y beneath the pond, i s the most important region to sample i n order to p r e d i c t d i s c h a r g e from the pond. Sampling i n region E, on the other hand, would seem to have q u e s t i o n a b l e v a l u e . For comparison, two other areas are COMPUTER SIMULATIONS / 55 0 32 64 96 128 160 192 224 256 288 320 HORIZONTAL DISTANCE, metres F i g . 8. L o c a t i o n of measurements i n regions E through H. a l s o t e s t e d ; measurements are made i n the c e n t r a l region (G) and near the outflow boundary (r e g i o n H). The next set of c o n d i t i o n a l s i m u l a t i o n s w i l l address the t h i r d q u e s t i o n . Once the l o c a t i o n f o r the i n i t i a l borehole has been determined, the s i t e i n v e s t i g a t o r may wish to c o l l e c t samples i n another l o c a t i o n . Four sampling s t r a t e g i e s are t e s t e d (schemes I, J , K, and L) that each i n c l u d e the i n i t i a l borehole l o c a t i o n chosen from the f i r s t set of c o n d i t i o n a l s i m u l a t i o n s and a second borehole (Figure 9) . F i g u r e 10 shows three p o s s i b l e sampling s t r a t e g i e s used •in a d d r e s s i n g the f o u r t h q u e s t i o n . In sampling s t r a t e g y M, s i x boreholes are d i s t r i b u t e d evenly throughout the flow COMPUTER SIMULATIONS / 56 F i g . 9. Sampling schemes I through L with two boreholes. COMPUTER SIMULATIONS / 57 F i g . 10. L o c a t i o n of measurements f o r sampling s t r a t e g i e s M through 0. domain with three measurements i n each borehole. Sampling s t r a t e g i e s N and 0 each have 30 measurements but d i f f e r i n COMPUTER SIMULATIONS / 58 the l o c a t i o n s . Both i n c l u d e the 18 measurements of s t r a t e g y M. In s t r a t e g y N, the measurements of s t r a t e g y M are augmented by four new boreholes with three measurements i n each. In s t r a t e g y 0, an a d d i t i o n a l two measurements are made in each of the s i x boreholes of s t r a t e g y M. These sampling s t r a t e g i e s are a p p l i e d to the s t r a t i f i e d r e f e r e n c e cases 14 through 17. Sampling schemes P through R, shown i n F i g u r e 11, address the l a s t set of q u e s t i o n s . In sampling scheme P, three measurements are taken i n each of four boreholes spaced evenly i n the flow domain. In scheme Q, three measurements are taken i n each of seven boreholes spaced evenly i n the domain. In scheme R, data are c o l l e c t e d i n 13 boreholes that have three measurements each. A summary of the t e s t runs i s presented i n Table 2. COMPUTER SIMULATIONS / 59 Scheme Q o -\ 1 1 1 1 1 1 1 1 1 O 32 64 96 128 160 192 224 256 288 320 HORIZONTAL DISTANCE, metres F i g . 11. L o c a t i o n of measurements s t r a t e g i e s P through R. f o r sampling COMPUTER SIMULATIONS / 60 TABLE 2. C o n d i t i o n a l S i m u l a t i o n s Reference case number Sampling strategy A - D E-H " I-L M - 0 P-R 1-5 using Sy X 1-5 using X X C-10 using S™ X X 11-13 using S™ X X 14-17 using S™ X Conditional simulation input Input to the c o n d i t i o n a l s i m u l a t i o n model f o r each sampling s t r a t e g y c o n s i s t s of (1) the number of measurements, the l o c a t i o n of measurements, the measured va l u e s of h y d r a u l i c c o n d u c t i v i t y as taken from the s p e c i f i c b l o c k s of the r e f e r e n c e case, (2) the standard d e v i a t i o n , and (3) the h o r i z o n t a l and v e r t i c a l i n t e g r a l s c a l e s of the h y d r a u l i c c o n d u c t i v i t y f i e l d . For each sampling s t r a t e g y a p p l i e d to a p a r t i c u l a r r e f e r e n c e case, a c o n d i t i o n a l s i m u l a t i o n comprising 300 runs i s generated and p r o b a b i l i t y d i s t r i b u t i o n s on estimates of pond d i s c h a r g e computed. Although the number of Monte C a r l o runs was s e l e c t e d a r b i t r a r i l y , Table 3 suggests that model output s t a b i l i z e s i n l e s s than 300 runs f o r both u n c o n d i t i o n a l and c o n d i t i o n a l COMPUTER SIMULATIONS / 61 TABLE 3. Comparison of Sim u l a t i o n s with 200, 300, and 400 Monte C a r l o R e a l i z a t i o n s Unconditional simulation (A, = 32m, A z = 12m) MG Input S" V Q 5, My* 200 -3.0 0.5 -3.00 0.48 0.65 0.51 0.78 0.23 0.26 300 -3.0 0.5 -3.00 0.43 0.63 0.48 0.76 0.23 0.23 400 -3.0 0.5 -3.00 0.48 0.63 0.4S 0.76 0.23 0.23 Conditional simulation using samples collected in region G for reference case 6 (/(y = -3.00, (Ty = 0.5, A, = 32m, Xz = 12m) A/C Input Yc Q Q-Qt RMSE (Syc)2 Y * 0 y 200 -2.74 0.32 -2.75 0.97 0.58 0.71 0.41 0.10 0.1S 300 -2.74 0.32 -2.75 0.03 0.57 0.70 0.42 0.10 0.16 400 -2.74 0.32 -2.75 0.06 0.57 0.70 0.42 0.10 0.16 * fs, Oy, Y, Sy for K in m/«; Q, Qt, Sg, RMSE in (m2/a x 10 - 2) s i m u l a t i o n s . In p a r t i c u l a r , the average v a r i a n c e i n l o g h y d r a u l i c c o n d u c t i v i t y over the s i m u l a t i o n , (S^)2, and the average v a r i a n c e i n estimates of pond d i s c h a r g e over the s i m u l a t i o n , ( S g ) 2 , are s i m i l a r for 200, 300, and 400 r e a l i z a t i o n s when a^=0.5. When the number of measurements i s s m a l l , i t i s d i f f i c u l t to a c c u r a t e l y estimate the standard d e v i a t i o n from the r e f e r e n c e case, Sr . For t h i s reason c o n d i t i o n a l s i m u l a t i o n s f o r the s t r a t e g i e s with 12 measurements (A, B, C, D) are generated two d i f f e r e n t ways. In one s i m u l a t i o n COMPUTER SIMULATIONS / 62 the standard d e v i a t i o n of the sample, S™, i s input to the c o n d i t i o n a l model and i n the other, Sr i s used. In t h i s way the i n f l u e n c e of the accuracy of estimates of the standard d e v i a t i o n (input parameter) on the accuracy of d i s c h a r g e estimates (output) i s examined. The number of data p o i n t s used in t h i s study i s i n s u f f i c i e n t to p r o p e r l y estimate e i t h e r the h o r i z o n t a l and v e r t i c a l i n t e g r a l s c a l e s or the form of the autocovariance f u n c t i o n . For t h i s reason X , X^ and the ex p o n e n t i a l form of (1) are assumed known. As mentioned above e i t h e r S™ or Sry i s used f o r a 2 i n ( 1 ) . These three parameters (X^, X z and standard d e v i a t i o n ) are used to set up the autocovariance matrix used i n forming the k r i g i n g e s t i m a t e , Y (i,j), and the 300 c o n d i t i o n a l r e a l i z a t i o n s , s(i,j). Conditional simulation output The output from each c o n d i t i o n a l s i m u l a t i o n c o n s i s t s of the average value and standard d e v i a t i o n i n log h y d r a u l i c c o n d u c t i v i t y formed over the set of 300 r e a l i z a t i o n s . These valu e s are denoted Y Q and Scy because they r e f e r to c o n d i t i o n a l s i m u l a t i o n s . In a d d i t i o n an estimate of pond d i s c h a r g e i s computed from the stream f u n c t i o n s o l u t i o n f o r each r e a l i z a t i o n . The average value of the estimate of pond d i s c h a r g e , Q, and standard d e v i a t i o n , S , are c a l c u l a t e d f o r each c o n d i t i o n a l s i m u l a t i o n . Combining these estimates with COMPUTER SIMULATIONS / 63 the a c t u a l d i s c h a r g e from the pond, Q{, and a co s t c o e f f i c i e n t , , the l o s s f u n c t i o n (13) can be c a l c u l a t e d . C O S T C O E F F I C I E N T S Table 4 shows the c o s t c o e f f i c i e n t s used i n c a l c u l a t i n g the sampling c o s t s f o r each s t r a t e g y . The sampling c o s t s are wi t h i n ten percent of those reported by Nakamoto et al [1986]. The c o s t c o e f f i c i e n t f o r the l o s s f u n c t i o n , , depends upon the cost of a l t e r n a t i v e treatment, pond width, and expected l i f e of the pond Cl ~ Lw Ct Th where C[ cost c o e f f i c i e n t [$T/L2]; Lw width of the pond [ L ] ; Ct c o s t of a l t e r n a t i v e treatment [$/L3]; Th time h o r i z o n [T], For a l l s i m u l a t i o n s , the le n g t h (76 m) and width (10 m) of the pond are f i x e d . A time h o r i z o n of 5 years i s used f o r a l l the s i m u l a t i o n s . When c a l c u l a t i n g the c o s t c o e f f i c i e n t , no account i s made f o r changes i n the value of the d o l l a r COMPUTER SIMULATIONS / 64 TABLE 4. Cost C o e f f i c i e n t s Symbol Description Cost (1986 $U.S.)* c d drilling 65.50/m drilled c e piezometer casing 1.30/m cased c v miscellaneous costs per saturated sample 508.20/piezometer c. miscellaneous costs per unsaturated sample 160/unsaturated sample c h miscellaneous costs per borehole 850/hole Cm miscellaneous costs per strategy 2,500 * Based on price lists and conversations with Westbay Instruments, Inc.; R. S. Technical Instruments, Ltd.; and Klohn Lconoff Consulting Engineers, all of Vancouver, Canada over the time h o r i z o n . Because i t i s d i f f i c u l t to determine the cost of a l t e r n a t i v e treatment (C ?) and the cost c o e f f i c i e n t (C{) f o r a h y p o t h e t i c a l case, the s e n s i t i v i t y of the r e s u l t s to d i f f e r e n t values w i l l be d i s c u s s e d . The design of sampling programs can be expected to vary depending upon the value of C|, a measure of the economic l o s s e s a s s o c i a t e d with poor p r e d i c t i o n s . Using an a r b i t r a r y value of , the value of the o b j e c t i v e f u n c t i o n , Z, i s c a l c u l a t e d , the sampling s t a t e g i e s compared, and the s e n s i t i v i t y d i s c u s s e d . COMPUTER TIME Computer s i m u l a t i o n s were c a r r i e d out on an FPS-164/MAX att a c h e d processor at the U n i v e r s i t y of B r i t i s h Columbia. For a c o n d i t i o n a l s i m u l a t i o n with 300 r e a l i z a t i o n s , the COMPUTER SIMULATIONS / 65 computer time was about 12 minutes. RESULTS HOMOGENEOUS CASE In order to gain a b e t t e r understanding of how h e t e r o g e n e i t i e s a f f e c t d i s c h a r g e from a pond, l e t us f i r s t examine the homogeneous case. F i g u r e 12a i s a contour p l o t of the stream f u n c t i o n s o l u t i o n f o r a uniform h y d r a u l i c c o n d u c t i v i t y . The contour l i n e s c o i n c i d e with s t r e a m l i n e s . T o t a l flow through each stream tube i s the same. Where the st r e a m l i n e s are c l o s e l y spaced, the h y d r a u l i c g r a d i e n t and groundwater v e l o c i t y are hi g h . The upper region of flow, t h e r e f o r e , i s more a c t i v e than the lower r e g i o n because the v e l o c i t i e s are hi g h e r . Note a l s o that the h y d r a u l i c g r a d i e n t beneath the pond i s nonuniform, reaching a maximum at the downstream edge. F i g u r e 12b i s a p l o t of f l u i d f l u x normal to the boundary for the case where K = 10" 3 m/s. Values above the h o r i z o n t a l l i n e represent recharge across a boundary segment and values below the l i n e represent d i s c h a r g e . Note that the recharge r a t e below the pond i s nonuniform and i n c r e a s e s to a maximum below the downstream edge of the pond. The higher v e l o c i t y and h y d r a u l i c g r a d i e n t i n t h i s area allow f o r gr e a t e r i n f i l t r a t i o n . Discharge from the pond i s 0.44 x 10~ 2 m 2/s. Between the pond and the outflow face, the f l u x across the water t a b l e i n the homogeneous case i s zero. 66 RESULTS / 67 F i g . 12. Homogeneous case with Y = -3. (a) Streamline p l o t with contour i n t e r v a l 8.0 x 10"' m2/s. (b) Flux p r o f i l e . The v e r t i c a l dashed l i n e in t h i s and subsequent f i g u r e s r e p r e s e n t s the i n t e r s e c t i o n between the water t a b l e and the v e r t i c a l boundary at the r i g h t hand s i d e . To the r i g h t of t h i s l i n e , the normal f l u x becomes h o r i z o n t a l d i s c h a r g e . RESULTS / 68 A PRIORI (UNCONDITIONAL) MODEL P r i o r to a n a l y z i n g the r e s u l t s of the c o n d i t i o n a l s i m u l a t i o n s , i t i s important to examine the s e n s i t i v i t y of the a p r i o r i heterogeneous model to changes i n the parameters that d e s c r i b e s p a t i a l v a r i a b i l i t y . Table 5 demonstrates the i n f l u e n c e of a , a , X , and X on y y x z d i s c h a r g e s t a t i s t i c s f o r u n c o n d i t i o n a l s i m u l a t i o n s . As expected, an order of magnitude decrease i n mean h y d r a u l i c c o n d u c t i v i t y causes an order of magnitude decrease in the average estimate of pond d i s c h a r g e , Q, and standard d e v i a t i o n , S^. The c o e f f i c i e n t of v a r i a t i o n , (u = / Q), a measure of the v a r i a b i l i t y of d i s c h a r g e e s t i m a t e s , remains the same. As shown on the middle segment of Table 5, the grea t e r the standard d e v i a t i o n i n h y d r a u l i c c o n d u c t i v i t y , a , the gr e a t e r the average value and standard d e v i a t i o n of pond d i s c h a r g e . One would expect the v a r i a b i l i t y i n discharge to in c r e a s e as the standard d e v i a t i o n of h y d r a u l i c c o n d u c t i v i t y i n c r e a s e s , but the f a c t o r s that e x p l a i n the i n c r e a s e i n average d i s c h a r g e are more complex. To understand t h i s behavior, c o n s i d e r the frequency d i s t r i b u t i o n s shown s c h e m a t i c a l l y i n F i g u r e 13. The frequency d i s t r i b u t i o n s f o r both K and l o g K are shown f o r a low standard d e v i a t i o n and fo r a high standard d e v i a t i o n . As the standard d e v i a t i o n in h y d r a u l i c c o n d u c t i v i t y i n c r e a s e s , the lognormal d i s t r i b u t i o n RESULTS / 69 TABLE 5. Example S e n s i t i v i t i e s f o r Discharge S t a t i s t i c s Based on the A P r i o r i Model Mean hydraulic conductivity Q* - 2 4.58 0.995 .22 - 3 0.457 0.105 .23 -4 0.046 0.010 .21 {(Ty = 0.2, Xx = 32, Xz = 12) Hydraulic conductivity standard deviation Q sq vi 0.2 0.457 0.105 .23 0.5 0.628 0.479 .76 0.9 1.33 2.72 2.0 0*5, = - 3 , Xx = T.2, X, = 12) Hydraulic conductivity correlation Aj.:A- Xx Xz Q 0 0 0.445 0.026 0.06 2.67 10 6 0.456 0.066 0.14 2.C7 32 12 0.457 0.105 0.23 2.67 48 18 0.464 0.124 0.27 2.67 16 6 0.456 0.066 0.14 5.33 32 6 0.461 0.079 0.17 8.00 48 6 0.470 0.090 0.19 5.33 32 6 0.461 0.079 0.17 2.67 32 12 0.457 0.105 0.23 1.78 32 18 0.451 0.114 0.25 (HV = - 3 , (uy + a ). Mid K : ( M " oy) < Y i < (M^ + a^). Low K : Y i < (n y - oy). (b) Streamline p l o t with contour i n t e r v a l 8.0 x 10-" m 2/s. (c) Flux p r o f i l e . RESULTS / 80 HORIZONTAL DISTANCE, metres C. Fluid Flux Normal To Surface ( 1 0 ~ 4 m/s) - i t i Edge of pond " J Y RESULTS / 81 flow of e f f l u e n t through the pond i s much l e s s than i f these zones were l o c a t e d d i r e c t l y beneath the pond. In f a c t , a high h y d r a u l i c c o n d u c t i v i t y zone i s l o c a t e d d i r e c t l y beneath the pond. The r e s u l t i s that pond d i s c h a r g e i n t h i s r e f e r e n c e case i s higher than i n the other four r e f e r e n c e cases generated from the same a p r i o r i model. F i g u r e 156 shows the s t r e a m l i n e s f o r r e f e r e n c e case 1 using a contour i n t e r v a l of 8.0 x 10"' m 2/s. As i n the homogeneous case, t o t a l flow through each stream tube i s the same. The region of most a c t i v e flow i s where the stream tubes are the densest and the region of l e a s t a c t i v e flow i s where the stream tubes are the most spread out. To a c e r t a i n extent, the regions of g e n e r a l l y a c t i v e or i n a c t i v e flow are set by the geometry of the flow domain. The lower region of flow tends to be i n a c t i v e and the upper r e g i o n of flow, p a r t i c u l a r l y i n the v i c i n i t y of the downstream edge of the pond, tends to be more a c t i v e . The d i f f e r e n c e s i n the s t r e a m l i n e s from the homogeneous case are caused by the h e t e r o g e n e i t i e s i n K. In g e n e r a l , the e f f e c t of low K zones i n impeding flow can be observed as regions where the s t r e a m l i n e s spread out. The small area of low K near the water t a b l e 50 m downstream of the pond causes water to e x i t a c r o s s the water t a b l e . The f l u x p r o f i l e f o r r e f e r e n c e case 1 i s shown i n F i g u r e 15c. The recharge area on the l e f t r e p r e s e n t s the RESULTS / 82 pond. As i n the homogeneous case, flow through the pond i n c r e a s e s from the l e f t edge of the pond at the impermeable boundary to a maximum near the downstream edge of the pond. Because of the high h y d r a u l i c c o n d u c t i v i t y zones u n d e r l y i n g the pond and the tendency of low h y d r a u l i c c o n d u c t i v i t y zones to be near the impermeable boundaries, flow through the pond i s higher i n r e f e r e n c e case 1 than i n the homogeneous case and the other four r e f e r e n c e cases. The di s c h a r g e area to the r i g h t of the v e r t i c a l dashed l i n e r e p r e s e n t s h o r i z o n t a l f l u x a c r o s s the r i g h t s i d e boundary. Flux i s nonuniform along t h i s boundary due to the he t e r o g e n e i t i e s . An important d i f f e r e n c e between the f l u x p r o f i l e s f o r the heterogeneous cases and the homogeneous case i s i n the p a t t e r n s of recharge and d i s c h a r g e a c r o s s the water t a b l e . Because a f r e e s u r f a c e model was used to f i n d the water t a b l e i n the homogeneous case, there should be no f l u x a c r o s s i t t h e o r e t i c a l l y . The small f l u x that does occur i s a r e s u l t of d i s c r e t i z a t i o n e r r o r . In the heterogeneous s i m u l a t i o n s , however, the water t a b l e i s f i x e d . Recharge a c r o s s the water t a b l e occurs i n areas where the water t a b l e would f a l l i f a f r e e s u r f a c e model was used. S i m i l a r l y , d i s c h a r g e across the water t a b l e r e p r e s e n t s areas where the water t a b l e would r i s e i f not h e l d f i x e d . I f a low K zone occurs near the water t a b l e , as i n r e f e r e n c e case 1, some of RESULTS / 83 the pond e f f l u e n t avoids t h i s zone by e x i t i n g a c r o s s the water t a b l e . T h i s may r e s u l t i n g r e a t e r pond di s c h a r g e than in a f r e e s u r f a c e model that a l l o w s the water t a b l e to r i s e and f o r c e s the e f f l u e n t through t h i s zone. It can be argued that c a l c u l a t i n g f l u x through the pond using a f i x e d water t a b l e might not be as r e a l i s t i c as c a l c u l a t i o n s using a f r e e s u r f a c e model. A s t o c h a s t i c f r e e s u r f a c e model, however, i s beyond the scope of the present study. P r o v i d e d the l i m i t a t i o n s of using a f i x e d water t a b l e are kept i n mind, i t i s a reasonable a l t e r n a t i v e to a s t o c h a s t i c f r e e s u r f a c e model. In Monte C a r l o s i m u l a t i o n s that use a f i x e d water t a b l e , a recharge area for one r e a l i z a t i o n may be a d i s c h a r g e area f o r another r e a l i z a t i o n . The water t a b l e from the homogeneous case was chosen as a reasonable average l o c a t i o n f o r the i n f i n i t e number of p o s s i b l e water t a b l e c o n f i g u r a t i o n s f o r the heterogeneous cases. The f l u x p r o f i l e f o r r e f e r e n c e case 1 r e f l e c t s the p a r t i c u l a r arrangement of h y d r a u l i c c o n d u c t i v i t y v a l u e s . A high K zone immediately downstream of the pond i s adjacent to a s m a l l low K zone near the water t a b l e . T h i s causes some of the pond e f f l u e n t to a v o i d the low K zone by e x i t i n g a c r o s s the water t a b l e on the upstream s i d e . Recharge occurs immediately downstream of t h i s zone. The d i s c h a r g e at the toe of the slope occurs where the high K zone i s near the RESULTS / 84 water t a b l e . F i g u r e 16 shows maps of l o g h y d r a u l i c c o n d u c t i v i t y , s t r e a m l i n e p l o t s , and f l u x p r o f i l e s f o r r e f e r e n c e cases 2 through 5. In r e f e r e n c e case 2, the pond i s u n d e r l a i n p r i m a r i l y by h y d r a u l i c c o n d u c t i v i t y values i n the mid-range. The high K v a l u e s appear to be more d i s p e r s e d than i n r e f e r e n c e case 1, with the e x c e p t i o n of a zone near the outflow boundary on the r i g h t hand s i d e . T h i s high K zone does not appear to have an e f f e c t on flow through the pond, probably because of i t s l o c a t i o n . As i n r e f e r e n c e case 1, a r e g i o n of low K i s near the lower l e f t corner i n a l e s s a c t i v e area of flow with g e n e r a l l y low v e l o c i t i e s . T h i s low K region has l i t t l e i n f l u e n c e on pond d i s c h a r g e . An area of mixed low and mid K occurs midway through the flow domain. Because most of the flow l i n e s from the pond pass through t h i s zone, d i s c h a r g e from the pond i s lower than i n r e f e r e n c e case 1. A small volume of water near the upper boundary can a v o i d t h i s zone, however, by e x i t i n g a c r o s s the water t a b l e , as seen i n the s t r e a m l i n e p l o t . The contour i n t e r v a l between s t r e a m l i n e s i s the same in t h i s f i g u r e as i n the p r e v i o u s f i g u r e . The r e s u l t of the s p a t i a l arrangement of K i n r e f e r e n c e case 2 i s that Q{ i s c l o s e r to the average d i s c h a r g e f o r the u n c o n d i t i o n a l s i m u l a t i o n than Qt f o r r e f e r e n c e case 1. In r e f e r e n c e case 3, there i s a h i g h h y d r a u l i c RESULTS / 85 F i g . 16. Reference cases 2 through 5. (a) Map of l o g h y d r a u l i c c o n d u c t i v i t y . (b) Streamline p l o t . Contour i n t e r v a l i s 8.0 x l0-«m 2/s. (c) Flux p r o f i l e . Fluid Flux Normal to Surface <10"4m/3) Reference Case 2 »on Reference Case 3 -8.0 J Reference Case 4 Reference Case 5 HORIZONTAL DISTANCE, metres _ CO RESULTS / 87 c o n d u c t i v i t y zone immediately downstream of the pond with a low h y d r a u l i c c o n d u c t i v i t y zone f u r t h e r downstream from i t . T h i s arrangement of c o n t r a s t i n g h y d r a u l i c c o n d u c t i v i t y zones causes l a r g e volumes of pond water to e x i t along the slope before reaching the v e r t i c a l outflow boundary. Because f l u i d near the water t a b l e can avoid t h i s low K zone, flow through the pond i s l i t t l e a f f e c t e d by i t . The average h y d r a u l i c c o n d u c t i v i t y of the upper and more a c t i v e region of flow i s most l i k e l y higher than u , r e s u l t i n g i n a Qt that i s s l i g h t l y higher than the expected value from the ensemble s t a t i s t i c s . Reference case 4 has a lower average h y d r a u l i c c o n d u c t i v i t y than the other four r e f e r e n c e cases. There are fewer high K blocks and more low K b l o c k s . The hig h K blocks have l i t t l e e f f e c t on di s c h a r g e from the pond because of t h e i r l o c a t i o n . They are e i t h e r near the impermeable boundary i n a l e s s a c t i v e region of flow or near the water t a b l e at the r i g h t hand s i d e , downstream of a low K zone. The low K zones occur throughout the a c t i v e region of flow. P a r t i c u l a r l y c r i t i c a l i n i n f l u e n c i n g flow i s a low K zone l o c a t e d below the downstream edge of the pond. R e c a l l from the homogeneous case that t h i s area had the highest v e l o c i t i e s and h y d r a u l i c g r a d i e n t r e s u l t i n g in the highest i n f i l t r a t i o n r a t e . In r e f e r e n c e case 4, the presence of a low K zone in t h i s area impedes the r a t e of i n f i l t r a t i o n RESULTS / 88 through a p o r t i o n of the pond r e s u l t i n g i n the lowest value of Qt . T h i s i s evident by comparing the f l u x from the pond on the f i v e f l u x p r o f i l e s . Comparing the s t r e a m l i n e p l o t s , the lower volume of flow through the pond i s evident from the s m a l l e r number of stream tubes o r i g i n a t i n g from the pond. In r e f e r e n c e case 5, although a p o r t i o n of the upstream end of the pond i s u n d e r l a i n by a low h y d r a u l i c c o n d u c t i v i t y zone, a hig h h y d r a u l i c c o n d u c t i v i t y l a y e r e x i s t s i n the upper three meters beneath the pond. The arrangement of high h y d r a u l i c c o n d u c t i v i t y zones i n the remainder of the flow domain r e s u l t s i n a l a r g e r volume of water fl o w i n g through the pond than the average value of d i s c h a r g e from the a p r i o r i model. Reference cases 6 through 16 have a l a r g e r standard d e v i a t i o n (oy = 0.5) and a v a r i e t y of i n t e g r a l s c a l e s . They w i l l be d i s c u s s e d i n l a t e r s e c t i o n s when the r e s u l t s p e r t a i n i n g to those r e f e r e n c e cases are d i s c u s s e d . In the next s e c t i o n , the technique of c o n d i t i o n a l s i m u l a t i o n w i l l be b r i e f l y reviewed while d e s c r i b i n g the r e s u l t s of the f i r s t s i m u l a t i o n s . RESULTS / 89 LOCATION FOR I N I T I A L BOREHOLES The f i r s t set of s i m u l a t i o n s i s designed to address the q u e s t i o n of where to l o c a t e i n i t i a l boreholes d u r i n g a s i t e i n v e s t i g a t i o n . The sampling s t r a t e g i e s d i s c u s s e d p r e v i o u s l y ( F i g u r e 7) are a p p l i e d to r e f e r e n c e cases 1 through 5. The r e s u l t s are presented in Table 8. The columns c o n t a i n summary s t a t i s t i c s on the porous medium, volume d i s c h a r g e , d e v i a t i o n s from Qf , the root mean square e r r o r and c o e f f i c i e n t of v a r i a t i o n . The rows are grouped a c c o r d i n g to borehole l o c a t i o n f o r each of the f i v e r e f e r e n c e cases. Estimates of the mean value and standard d e v i a t i o n of lo g K are obtained from the set of borehole measurements (Y m and S™). Because of s t a t i s t i c a l f l u c t u a t i o n s due to the number and l o c a t i o n of measured v a l u e s , Ym and S™ w i l l d i f f e r from the average value and standard d e v i a t i o n of l o g K from the r e f e r e n c e case, Yr and Sry. In order to i n v e s t i g a t e the e f f e c t of borehole l o c a t i o n and not the e f f e c t of e r r o r s i n e s t i m a t i n g the standard d e v i a t i o n of l o g K, the standard d e v i a t i o n of the r e f e r e n c e case, Sry, was used as input to the c o n d i t i o n a l model r a t h e r than the standard d e v i a t i o n of the measurements, S™. L a t e r i n t h i s s e c t i o n , r e s u l t s w i l l be presented from s i m u l a t i o n s using V The estimate, 7 . and standard d e v i a t i o n are used as input to the c o n d i t i o n a l model (columns 1 and 2, Table 8), RESULTS / 90 TABLE 8. C o n d i t i o n a l S i m u l a t i o n s of Reference Cases 1 through 5 using S i n g l e Boreholes and Assuming the Standard D e v i a t i o n of the Reference Case i s Known RCN y * Qt Q Q-Qt 5< RMSE Borehole A 1 -2.84 0.20 -2.89 0.20 0.73 0.62 -0.11 0.15 0.19 0.25 2 -2.94 0.17 -2.96 0.16 0.46 0.50 0.04 0.09 0.10 0.19 3 -3.02 0.20 -3.02 0.19 0.48 0.44 -0.04 0.11 0.11 0.24 4 -3.08 0.18 -3.08 0.17 0.33 0.38 0.05 0.07 0.09 0.19 5 -3.16 0.19 -3.08 0.19 0.51 0.39 -0.13 0.09 0.15 0.23 Average 0.19 0.18 0.07** 0.10 0.13 Std Dev 0.01 0.01 0.04 0.03 0.04 Borehole B 1 -2.71 0.20 -2.71 0.19 0.73 0.95 0.22 0.16 0.27 0.17 2 -2.90 0.17 -2.90 0.17 0.46 0.54 0.08 0.09 0.12 0.17 3 -2.98 0.20 -3.04 0.19 0.48 0.42 -0.05 0.08 0.10 0.19 4 -3.30 0.18 -3.27 0.17 0.33 0.23 -0.11 0.04 0.11 0.16 5 -2.99 0.19 -2.97 0.18 0.51 0.51 -0.00 0.08 0.08 0.16 Average 0.19 0.18 0.09 0.09 0.14 Std Dev 0.01 0.01 0.08 0.05 0.08 Borehole C 1 -2.74 0.20 -2.77 0.19 0.73 0.77 0.04 0.20 0.21 0.26 2 -2.97 0.17 -2.96 0.16 0.46 0.51 0.05 0.10 0.11 0.20 3 -2.81 0.20 -2.80 0.19 0.48 0.73 0.26 0.18 0.31 0.24 4 -3.14 0.18 -3.20 0.17 0.33 0.29 -0.05 0.06 0.08 0.21 5 -2.92 0.19 -2.83 0.19 0.51 0.69 0.17 0.15 0.23 0.22 Average 0.19 0.18 0.11 0.14 0.19 Std Dev 0.01 0.01 0.10 0.06 0.09 Borehole D 1 -2.82 0.20 -2.87 0.19 0.73 0.61 -0.12 0.16 0.20 0.26 2 -2.74 0.17 -2.80 0.17 0.46 0.72 0.26 0.14 0.30 0.19 3 -3.06 0.20 -3.06 0.19 0.48 0.42 -0.06 0.11 0.12 0.26 4 -2.96 0.18 -2.05 0.17 0.33 0.50 0.16 0.10 0.19 0.20 5 -2.94 0.19 -2.91 0.19 0.51 0.60 0.08 0.15 0.17 0.24 Average 0.19 0.18 0.14 0.13 0.20 Std De\ 0.01 0.01 0.08 0.02 0.06 * Y, S„ for K in m/a; Q, RMSE in ( m 2 / * x 10~2) **E[\Q--Qt\\ RESULTS / 91 along with the set of measured K v a l u e s . The output from a c o n d i t i o n a l s i m u l a t i o n i n c l u d e s the average value and standard d e v i a t i o n of l o g K {Y c and Scy) formed over the set of r e a l i z a t i o n s (columns 3 and 4). Because of the f i n i t e number of r e a l i z a t i o n s used i n the a n a l y s i s , the input (Ym and standard d e v i a t i o n ) w i l l d i f f e r from the output (Y and Sc) A c o n d i t i o n a l Monte C a r l o s i m u l a t i o n i s c a r r i e d out by p r e s e r v i n g the values of h y d r a u l i c c o n d u c t i v i t y at the measurement l o c a t i o n s and g e n e r a t i n g 300 r e a l i z a t i o n s . For each r e a l i z a t i o n an estimate of d i s c h a r g e from the pond, Q, i s c a l c u l a t e d from the stream f u n c t i o n s o l u t i o n . Each r e a l i z a t i o n w i l l have a d i f f e r e n t Q because the random p a t t e r n s of s p a t i a l v a r i a b i l i t y c r e a t e unique flow p a t t e r n s . The average value of pond di s c h a r g e f o r the s i m u l a t i o n , Q, and the standard d e v i a t i o n , S q l are c a l c u l a t e d (columns 6 and 8, Table 8 ) . Sq i n d i c a t e s the v a r i a b i l i t y i n Q. Q i s the most l i k e l y value of pond d i s c h a r g e . Bias r e f e r s to the d i f f e r e n c e between the average value of an e s t i m a t o r and the parameter i t i s supposed to estimate [Benjamin and Cornell, 1970; Davis, 1973]. The b i a s i n the estimate (column 7) i s the d i f f e r e n c e between Q and the a c t u a l flow through the pond from the r e f e r e n c e case, Qt (column 5). The root mean square e r r o r (column 9) p r o v i d e s a b e t t e r measure of the u n c e r t a i n t y of p r e d i c t i o n s of flow through the pond than the RESULTS / 92 b i a s or the c o e f f i c i e n t of v a r i a t i o n because i t accounts f o r both the b i a s and the v a r i a b i l i t y i n pond di s c h a r g e (equation 12). The smaller the root mean square e r r o r (RMSE), the b e t t e r the p r e d i c t i o n . Using t h i s framework, the r e d u c t i o n i n p r e d i c t i o n u n c e r t a i n t y f o r each sampling s t r a t e g y , as measured by the RMSE, i s compared. By r e p e a t i n g the above steps f o r d i f f e r e n t r e f e r e n c e cases, the s e n s i t i v i t y of the r e s u l t s to d i f f e r e n c e s in the p a t t e r n s of s p a t i a l v a r i a b i l i t y can be examined. Because the o b j e c t i v e f u n c t i o n depends s t r o n g l y upon the values chosen for the cost c o e f f i c i e n t s , a d i s c u s s i o n of the economic i s s u e s i s postponed u n t i l l a t e r . Sy used as input Consider the r e s u l t s presented in Table 8 f o r r e f e r e n c e case 1 u s i n g the standard d e v i a t i o n of the r e f e r e n c e case, Sry, as i n p u t . A c o n d i t i o n a l s i m u l a t i o n based on K measurements from borehole C l o c a t e d midway between the pond and the r i g h t s i d e boundary y i e l d s a value of Q (0.77 x 10~ 2 m 2/s) t h a t i s c l o s e r to Qt (0.73 x 10" 2 m 2/s) than the s i m u l a t i o n s based on measurements below the pond (boreholes A or B) or near the r i g h t s i d e boundary (D). The average value of K i n borehole B i s higher than the other Ym and y i e l d s the highest Q. The two lowest Ym are from boreholes A and D, y i e l d i n g the lowest average estimates of d i s c h a r g e . RESULTS / 93 R e c a l l t h at the arrangement of h y d r a u l i c c o n d u c t i v i t y i n t h i s r e f e r e n c e case i s such that the upper region of flow c o n s i s t s of high and moderate K v a l u e s , while the low K v a l u e s are mostly c o n f i n e d to l e s s a c t i v e regions of flow. Qt i s q u i t e high i n t h i s r e f e r e n c e case and a f a i r l y high value of Ym i s needed to o b t a i n a small b i a s . The v a r i a b i l i t y i n estimates of pond d i s c h a r g e (S = 0.20) i s higher f o r the s i m u l a t i o n using borehole C than f o r the other s i m u l a t i o n s . The c o n d i t i o n a l s i m u l a t i o n based on measurements below the pond (A) r e s u l t s i n an underestimate of pond di s c h a r g e {Q = 0.62 x 10" 2 m 2/s), but the lowest v a r i a b i l i t y (S = 0.15) of the four s i m u l a t i o n s . The r e s u l t i s a root mean square e r r o r , which accounts f o r both the b i a s and the v a r i a b i l i t y of the estimates, that i s lowest f o r measurements taken below the middle of the pond (RMSE = 0.19 x 10" 2 m 2/s) f o r r e f e r e n c e case 1. In r e f e r e n c e case 2, the b i a s i n estimates of discharge i s lowest f o r the s i m u l a t i o n using borehole A l o c a t e d underneath the middle of the pond, and the v a r i a b i l i t y i s lowest f o r the s i m u l a t i o n based on measurements from the borehole near the downstream edge of the pond (B). Again, the root mean square e r r o r i s lowest f o r the s i m u l a t i o n u s i n g K data c o l l e c t e d from the borehole that i s l o c a t e d below the middle of the pond (A). In ref e r e n c e cases 3 and 5, the root mean square e r r o r i s lowest f o r s i m u l a t i o n s RESULTS / 94 based on measurements made i n borehole B (edge of pond). In re f e r e n c e case 4, although i s lowest f o r the s i m u l a t i o n based on measurements i n borehole B, the d e v i a t i o n of Q from Qt i s lowest f o r the s i m u l a t i o n using data from borehole C. With a low S^, the s i m u l a t i o n based on data from borehole C has the lowest root mean square e r r o r of the four s i m u l a t i o n s f o r r e f e r e n c e case 4. In add r e s s i n g the q u e s t i o n of where to l o c a t e i n i t i a l b o r eholes, we must keep in mind that the purpose of network design i s to do the best job, on average, of minimizing e r r o r . The best s t r a t e g y from a set of f i n i t e a l t e r n a t i v e s i s the s t r a t e g y that i s most l i k e l y to perform b e t t e r than the o t h e r s . Because of the v a r i a t i o n s i n s p a t i a l arrangement of h y d r a u l i c c o n d u c t i v i t y between re f e r e n c e cases, a best s t r a t e g y w i l l not be best f o r a l l r e f e r e n c e cases. Because an a c t u a l f i e l d s i t e i s analogous to a s i n g l e r e a l i z a t i o n or re f e r e n c e case, the r e s u l t s imply that a sampling s t r a t e g y that performs w e l l f o r a number of d i f f e r e n t s i t e s may perform q u i t e p o o r l y at a p a r t i c u l a r s i t e . In order to f i n d the p r o b a b i l i t y d i s t r i b u t i o n of the root mean square e r r o r f o r each borehole l o c a t i o n , a l a r g e number of s i m u l a t i o n s using m u l t i p l e r e f e r e n c e cases are r e q u i r e d . Although t h i s task i s not f e a s i b l e f o r the present study, we can make some o b s e r v a t i o n s from the r e s u l t s of the f i v e r e f e r e n c e cases that are co n s i d e r e d . In 2 of the 5 RESULTS / 95 re f e r e n c e cases, s i m u l a t i o n s using data c o l l e c t e d i n borehole A (below the middle of the pond) was the best s t r a t e g y of the four s t r a t e g i e s c o n s i d e r e d because they y i e l d e d the lowest RMSE. In 2 of the 5 r e f e r e n c e cases, borehole B ( l o c a t e d below the downstream edge of the pond) provided the lowest root mean square e r r o r . In r e f e r e n c e case 4, the root mean square e r r o r was lowest f o r s i m u l a t i o n s based on measurements i n borehole C. If the purpose of network design i s to minimize l e v e l 2 u n c e r t a i n t y , as measured by root mean square e r r o r , then the data suggest that l o c a t i n g i n i t i a l boreholes below the pond may be a b e t t e r s t r a t e g y than l o c a t i n g i n i t i a l boreholes elsewhere. These r e s u l t s apply to t h i s p a r t i c u l a r boundary value problem with l a r g e i n t e g r a l s c a l e s and small standard d e v i a t i o n , when Sr i s assumed known. By averaging the r e s u l t s over the f i v e r e f e r e n c e cases fo r each borehole, we can s p e c u l a t e f u r t h e r on the best l o c a t i o n f o r i n i t i a l b o r e h o l e s . The r e s u l t s are presented i n Table 8. Note that the average value and standard d e v i a t i o n of \Q - Qt |, not (Q - Qt ), i s presented because we are i n t e r e s t e d i n the magnitude of the d e v i a t i o n from true pond d i s c h a r g e , not whether the average Q underestimates or overestimates Qt . The average d e v i a t i o n i s lowest f o r s i m u l a t i o n s using data from borehole A, f o l l o w e d by s i m u l a t i o n s using data from boreholes B, C, and then D. RESULTS / 96 Although the average value and standard d e v i a t i o n of Scy i s n e a r l y the same f o r a l l boreholes, the average standard d e v i a t i o n i n dis c h a r g e estimates i s lowest from s i m u l a t i o n s u s i n g data c o l l e c t e d below the edge of the pond (B) and next lowest from s i m u l a t i o n s u s i n g data from below the middle of the pond (A). The average RMSE i s lowest from s i m u l a t i o n s using borehole A data, the next lowest average i s from s i m u l a t i o n s using borehole B data. The average RMSE i s lowest f o r s i m u l a t i o n s using data i n borehole C or borehole D. Of course, these comparisons must be viewed with c a u t i o n because of the small sample s i z e . S™ used as input The above s i m u l a t i o n s n e g l e c t the e r r o r s i n e s t i m a t i n g the standard d e v i a t i o n of the h y d r a u l i c c o n d u c t i v i t y d i s t r i b u t i o n , Sry. In order to i n c o r p o r a t e these e r r o r s , a l l of the t r i a l s d i s c u s s e d i n the pre v i o u s s e c t i o n that used Sry were repeated using S™ i n s t e a d (Table 9 ) . The d i f f e r e n c e between estimates of pond d i s c h a r g e when Sry i s assumed known and when i t i s estimated from the 12 measurements can be seen by comparing Tables 8 and 9. The average estimate of pond d i s c h a r g e , Q, i s i n g e n e r a l , l e s s s e n s i t i v e than S and RMSE to the value of standard d e v i a t i o n used as in p u t . Note that the estimates of the standard d e v i a t i o n i n h y d r a u l i c c o n d u c t i v i t y from the measurements, S™, tend to be lower RESULTS / 97 TABLE 9. C o n d i t i o n a l S i m u l a t i o n s of Reference Cases 1 through 5 using S i n g l e Boreholes and E s t i m a t i n g the Standard D e v i a t i o n from the Measurements RCN Y * I m cm Ye Qt Q Q-Qt sq RMSE Vq Borehole A 1 -2.84 0.17 -2.89 0.17 0.73 0.61 -0.12 0.13 0.i7 0.20 2 -2.94 0.07 -2.97 0.06 0.46 0.47 0.01 0.04 0.04 0.08 3 -3.02 0.10 -3.01 0.10 0.48 0.43 -0.04 0.05 0.07 0.12 4 -3.08 0.09 -3.10 0.09 0.33 0.38 0.04 0.03 0.06 0.09 5 -3.16 0.24 -3.08 0.23 0.51 0.40 -0.11 0.12 0.16 0.29 Average -3.01 0.13 0.07** 0.07 0.10 0.16 Standard Deviation 0.09 0.07 0.05 0.04 0.06 0.09 Borehole B 1 -2.71 0.16 -2.70 0.16 0.73 0.94 0.21 0.13 0.25 0.14 2 -2.90 0.25 -2.90 0.24 0.46 0.55 0.09 0.14 0.16 0.25 3 -2.98 0.16 -3.04 0.15 0.48 0.42 -0.06 0.06 0.09 0.15 4 -3.30 0.18 -3.26 0.12 0.33 0.23 -0.11 0.02 0.11 0.11 5 -2.99 0.15 -2.97 0.14 0.51 0.51 -0.01 0.06 0.06 0.12 Average -2.97 0.16 0.10 0.08 0.13 0.15 Standard Deviation 0.20 0.05 0.07 0.05 0.07 0.06 Borehole C 1 -2.74 0.18 -2.77 0.17 0.73 0.76 0.03 0.18 0.18 0.23 2 -2.97 0.12 -2.95 0.12 0.46 0.50 0.05 0.07 0.08 0.14 3 -2.81 0.14 -2.80 0.14 0.48 0.71 0.23 0.12 0.26 0.17 4 -3.14 0.14 -3.20 0.14 0.33 0.2S -0.05 0.05 0.07 0.16 5 -2.92 0.22 -2.83 0.21 0.51 0.70 0.18 0.18 0.25 0.25 Average -2.91 0.16 0.11 0.12 0.17 0.19 Standard Deviation 0.18 0.04 0.09 0.06 0.09 0.05 Borehole D 1 -2.82 0.14 -2.87 0.14 0.73 0.59 -0.14 0.10 0.18 0.18 2 -2.74 0.15 -2.80 0.15 0.46 0.72 0.26 0.12 0.28 0.17 3 -3.06 0.11 -3.06 0.11 0.48 0.40 -0.08 0.06 0.10 0.14 4 -2.96 0.17 -2.95 0.16 0.33 0.49 0.16 0.09 0.18 0.19 5 -2.04 0.22 -2.91 0.21 0.51 0.61 0.09 0.17 0.19 0.28 Average -2.92 0.15 0.15 0.11 0.19 0.19 Standard Deviation 0.10 0.04 0.07 0.04 0.07 0.05 * Y, S„ for A' in m/s; UJ — i LU 30 -l i U H I i l f <7y' < 0.85 S y i -< > LU 64 96 128 160 192 224 256 HORIZONTAL DISTANCE, metres 288 320 F i g . 18. Reduction of u n c e r t a i n t y i n h y d r a u l i c c o n d u c t i v i t y from sampling (a) i n region F and (b) i n borehole A. borehole A underestimate Yr and overestimate Sry while samples from region F overestimate Yr and underestimate Sry. The magnitude of the b i a s i s s i m i l a r f o r the two s i m u l a t i o n s but the lower S i n s i m u l a t i o n F leads to a lower RMSE. To summarize the r e s u l t s f o r refer e n c e cases 1 through 5, s i m u l a t i o n s c o n d i t i o n e d on 12 measurements i n a borehole below the middle of the pond l e a d , on the average, to lower RESULTS / 109 values of RMSE than s i m u l a t i o n s using 34 samples c o l l e c t e d in the upper 9 m below the pond. Sampling i n the d e p o s i t s immediately below the pond leads to lower average RMSE than sampling i n regions E, G, or H. These r e s u l t s apply to re f e r e n c e cases with small standard d e v i a t i o n i n K and l a r g e i n t e g r a l s c a l e s . In the next s e c t i o n s , sampling schemes E through H w i l l be a p p l i e d to r e f e r e n c e cases with l a r g e r standard d e v i a t i o n or smaller i n t e g r a l s c a l e s . Reference cases with l a r g e r standard d e v i a t i o n Reference cases 6 through 10 are generated using the l a r g e i n t e g r a l s c a l e s of the pr e v i o u s r e f e r e n c e cases and gr e a t e r standard d e v i a t i o n i n l o g K (o =0.5). The h y d r a u l i c c o n d u c t i v i t y maps and st r e a m l i n e p l o t s are presented i n F i g u r e 19. The f l u x p l o t s are presented i n F i g u r e 20. The contour i n t e r v a l i s the same f o r the f i v e s t r e a m l i n e p l o t s (2.0 x 10~ 3 m 2/s). I n c r e a s i n g the standard d e v i a t i o n c r e a t e s a l a r g e r c o n t r a s t between low K and high K zones. The r e s u l t i s more i r r e g u l a r flow p a t t e r n s with s t r e a m l i n e s spreading out to bypass low K zones and tending to conc e n t r a t e i n the high K zones. In r e f e r e n c e case 6, the pond i s u n d e r l a i n by low and mid K regions with a high K region at depth. Immediately downstream from the pond i s a l a r g e low K r e g i o n that i n t e r s e c t s the ground s u r f a c e and extends downward 40 m. RESULTS / 110 Reference Case 9 Reference Case 10 HORIZONTAL DISTANCE, metres HORIZONTAL DISTANCE, metres F i g . 19. Reference cases 6 through 10. (a) Log h y d r a u l i c c o n d u c t i v i t y maps, (b) Streamline p l o t s with contour i n t e r v a l 2.0 x 10~ 3 m 2/s. Fluid Flux Normal to Surface (10~4m/a) Reference Case 8 -10.0 10.0-1 Reference Case 7 o.o -5.0-RESULTS / 1 1 1 Reference Case 8 lo.o-i Reference Case 9 -5.0 -5.0--10 .0 J Reference Case 10 -i_r»' F i g . 20. Flux p r o f i l e s f o r ref e r e n c e cases 6 through 10, RESULTS / 112 T h i s low K region impedes the flow of water from the pond and the high K l a y e r at depth causes a l a r g e percentage of pond e f f l u e n t to t r a v e l deep i n the flow system. T h i s arrangement f o r c e s almost a l l of the pond e f f l u e n t to d i s c h a r g e i n t o the r i v e r and r e s u l t s i n the lowest value of Qt f o r r e f e r e n c e cases 6 through 10. T h i s can be seen by comparing the s t r e a m l i n e p l o t from r e f e r e n c e case 6 with the s t r e a m l i n e p l o t s from the other r e f e r e n c e cases. The number of stream tubes o r i g i n a t i n g below the pond i s lowest f o r r e f e r e n c e case 6, even though t o t a l flow through the system i s lowest f o r r e f e r e n c e case 10. The h i g h e s t value of Qt i s f o r r e f e r e n c e case 7 (Figure 21). E x t e n s i v e areas of high K e x i s t , with a l a r g e region extending from the pond deep i n t o the flow system. The high and mid K zones near the water t a b l e cause over 24% of the pond e f f l u e n t to e x i t a c r o s s the water t a b l e . Note that when Oy=0.5, the c o n t r a s t between high and low K zones r e s u l t s in more i r r e g u l a r flow p a t t e r n s . For example, i n the high K zone under the pond, a low K block i s surrounded by mid K b l o c k s . T h i s lower K area causes i n f i l t r a t i o n r a t e s d i r e c t l y above i t to be low. T h i s i s shown as a d i p i n the f l u x p l o t , that otherwise i n c r e a s e s from the l e f t hand boundary to the downstream edge of the pond. The high K l a y e r s that l i e along the b a s a l boundary in r e f e r e n c e case 8 draw much of the e f f l u e n t deep i n t o the RESULTS / 113 120 -\ 100 -co i ^ r ^ c o w ^ c o c o c \ i ^ o c » o ) c O i s - c o O C 0 « > a > C \ l i n C 0 i - ^ - r * . O C \ J U ) C 0 i - ^ > Q ( 1 0 ~ 2 m 2 / s ) F i g . 21. Frequency d i s t r i b u t i o n of pond d i s c h a r g e f o r the a p r i o r i model (vy=-3, ay=0.5, \x=32m, and .X =l2m). Pond d i s c h a r g e f o r r e f e r e n c e cases 6 through 10 shown f o r comparison. RESULTS / 1 1 4 flow system. Combined with the tendency of low K blocks to be i n the upper part of the flow domain, at l e a s t 99% of the pond e f f l u e n t i s f o r c e d to dis c h a r g e i n t o the r i v e r . The value of Qt i s between those of ref e r e n c e cases 6 and 7 and i s below average f o r the ensemble. Reference cases 9 and 10 behave q u i t e s i m i l a r l y . The pond i s u n d e r l a i n by low to moderate K zones. An e x t e n s i v e low K zone e x i s t s in the downstream p o r t i o n of the flow domain causing a l a r g e percentage of the pond e f f l u e n t to e x i t a c r o s s the water t a b l e before reaching the r i v e r . The high K region at depth does not appear to i n f l u e n c e d i s c h a r g e from the pond. Qt f o r these r e f e r e n c e cases i s lower than the expected value (0.63 x 10" 2 m 2/s) f o r the ensemble. S i m u l a t i o n s c o n d i t i o n e d on measurements i n regions E through H are generated f o r re f e r e n c e cases 6 through 10. The r e s u l t s are presented i n Table 12. Note that the behavior of i n d i v i d u a l r e a l i z a t i o n s i s much more v a r i a b l e than f o r the previous s i m u l a t i o n s with o^=0.2. The RMSE i s lowest f o r samples from r e g i o n E i n r e f e r e n c e case 6, f o r samples from region F i n r e f e r e n c e cases 7, 9, and 10, and fo r samples from region H i n re f e r e n c e case 8. The average value of RMSE i s lowest f o r samples from r e g i o n F. To understand t h i s behavior, l e t us examine a few cases i n more d e t a i l . RESULTS / 115 TABLE 12. C o n d i t i o n a l S i m u l a t i o n s of Measurements i n Regions E through H f o r Reference Cases 6 through 10 R C N Y * m ''y Qt Q Q-Qt S, RMSE Region E 6 -3.34 0.34 -3.27 0.33 0.38 0.28 -0.10 0.12 0.16 0.42 7 -3.29 0.35 -3.44 0.33 1.45 0.19 -1.26 0.09 1.26 0.46 8 -3.14 0.38 -3.02 0.37 0.55 0.49 -0.06 0.27 0.27 0.55 9 -3.C6 0.35 -3.51 0.33 0.41 0.17 -0.25 0.07 0.26 0.40 10 -3.22 0.40 -3.15 0.38 0.42 0.37 -0.05 0.21 0.21 0.55 Average 0.34** 0.15 0.43 Standard Deviation 0.52 0.09 0.46 Region F 6 -3.30 0.35 -3.54 0.34 0.38 0.19 -0.19 0.05 0.19 0.24 7 -2.65 0.39 -2.59 0.37 1.45 1.24 -0.21 0.34 0.40 0.27 8 -3.13 0.47 -3.18 0.46 0.55 0.29 -0.25 0.14 0.29 0.46 9 -3.21 0.25 -3.19 0.24 0.41 0.27 -0.14 0.06 0.15 0.20 10 -3.18 0.27 -3.07 0.26 0.42 0.35 -0.08 0.08 0.11 0.22 Average 0.17 0.13 0.23 Standard Deviation 0.07 0.12 0.11 Region G 6 -2.74 0.32 -2.75 0.31 0.38 0.95 0.57 0.40 0.70 0.42 7 -2.81 0.33 -2.75 0.32 1.45 0.94 -0.50 0.12 0.66 0.45 8 -3.36 0.30 -3.32 0.29 0.55 0.24 -0.31 0.09 0.32 0.39 9 -3.11 0.38 -3.08 0.38 0.41 0.50 0.09 0.24 0.25 0.47 10 -3.04 0.41 -3.04 0.30 0.42 0.64 0.22 0.24 0.33 0.38 Average 0.34 0.28 0.45 Standard Deviation 0.20 0.14 0.21 Region H 6 -2.95 0.29 -2.73 0.29 0.38 0.93 0.55 0.36 0.65 0.38 7 -2.80 0.40 -2.91 0.39 1.45 0.67 -0.78 0.35 0.86 0.53 8 -3.09 0.25 -3.05 0.24 0.55 0.41 -0.14 0.13 0.19 0.32 1 9 -2.88 0.26 -2.91 0.25 0.41 0.59 0.18 0.18 0.25 0.31 10 -3.29 0.38 -3.15 0.36 0.42 0.40 -0.02 0.22 0.22 0.55 Average 0.33 0.25 0.44 Standard Deviation 0.32 0.10 0.30 * Y, Sv for K in m/s; Q, RMSE in (m 2 / * x 1 0 - 2 ) **E[\Q-Qt\\ RESULTS / 1 1 6 Consider r e f e r e n c e case 6 . Because Qt i s c o n s i d e r a b l y lower than the average estimate of pond d i s c h a r g e from the u n c o n d i t i o n a l case, measurements that underestimate Yr (sample l o c a t i o n s E and F) are l i k e l y to have a lower b i a s than measurements that overestimate Y' (sample l o c a t i o n s G and H). I t i s i n t e r e s t i n g to note that i n t h i s p a r t i c u l a r case, the highest value of S™ (F) leads to the lowest value of S and the two lowest values of S™ (G and H) l e a d to the H y h i g h e s t values of 5^. Other f a c t o r s are more important, t h e r e f o r e , i n determining the v a r i a b i l i t y of est i m a t e s of pond d i s c h a r g e than the estimate of the standard d e v i a t i o n i n h y d r a u l i c c o n d u c t i v i t y . P r o p e r l y e s t i m a t i n g the standard d e v i a t i o n i n K should not be a prime o b j e c t i v e of a sampling program that i s designed to p r e d i c t d i s c h a r g e . In r e f e r e n c e case 7, the only low K zone i n a g e n e r a l l y h i g h K flow system i s l o c a t e d i n region E. Samples from here l e a d to a gross underestimate of Qt and a b i a s that i s h a l f an order of magnitude g r e a t e r than the b i a s from other s i m u l a t i o n s and a RMSE that i s h a l f an order of magnitude g r e a t e r . Because of the very high Qt , a l l the s i m u l a t i o n s underestimate Q(. The sampling s t r a t e g y with the hi g h e s t Ym (F) leads to the lowest value of RMSE. In summary, the behavior of the r e a l i z a t i o n s i s much more v a r i a b l e f o r the s i m u l a t i o n s based on these r e f e r e n c e cases than those based on the f i r s t f i v e r e f e r e n c e cases RESULTS / 1 1 7 because o y i s h i g h e r . For three of the f i v e r e f e r e n c e cases, sampling i n the region below the pond p r o v i d e d the lowest value of RMSE. The average value of RMSE was a l s o lowest f o r s i m u l a t i o n s based on samples c o l l e c t e d in t h i s r e g i o n . Reference cases with s m a l l e r i n t e g r a l s c a l e s S i m u l a t i o n s with measurements of K i n r e g i o n s E through H were repeated f o r r e f e r e n c e cases 11 through 13. These r e f e r e n c e cases have i n t e g r a l s c a l e s that are h a l f the s i z e of the p r e v i o u s i n t e g r a l s c a l e s and ay = 0.5. F i g u r e 22 shows the l o g K maps, s t r e a m l i n e p l o t s , and f l u x p r o f i l e s f o r the r e f e r e n c e cases. The contour i n t e r v a l i s the same fo r the three s t r e a m l i n e p l o t s (2.0 x 10~ 3 m 2/s). The pond i s u n d e r l a i n by a low K l a y e r at depth i n r e f e r e n c e case 11. Immediately downstream of the pond i s an e x t e n s i v e high K r e g i o n . As a r e s u l t , much of the pond e f f l u e n t remains in the upper p o r t i o n of the flow domain and at l e a s t 35% e x i t s a c r o s s the water t a b l e before r e a c h i n g the r i v e r . Discharge from the .pond does not seem to be a f f e c t e d by the low K regions at depth and i s the highest of the three r e f e r e n c e cases. In r e f e r e n c e case 12, the m a j o r i t y of the high K blocks are i n the upper l a y e r s and the m a j o r i t y of the low K blocks are i n the lower l a y e r s . Midway between the edge of the pond and the r i g h t s i d e boundary, a low K region extends to the ELEVATION, tn ELEVATION, m ELEVATION, m Ti to 3 T> rr » K> »-"< n> o Cb in n rt> • Ul 0> c re r - 1 3 « M • O it o (D tr >—• o c: 0) o n pu as TJ o c fD o o tn o nt ti H* • o < — » 1—' c — * fD t-t rt •< • • in rt rr 3 :r rt> 0> if TJ O < c OJ • lO 1—' o cr -M l w O J ro • o v i rr — i-l fi) X ID ^ B> 3 —• i—• O M . f • D O 81 I RESULTS / 1 1 9 water t a b l e . At l e a s t 63% of the e f f l u e n t from the pond flows through the shallow high K region and e x i t s across the water t a b l e upstream of the low K zone. The e f f e c t i s a high value of Qt (Figure 23). The lowest Qt i s from r e f e r e n c e case 13. The pond i s u n d e r l a i n by mid and low K regions with a high K region f u r t h e r below. The down s l o p i n g high K r e g i o n midway between the pond and the r i v e r draws the e f f l u e n t down i n t o the lower p o r t i o n of the flow domain causing at l e a s t 86% of the e f f l u e n t to discharge i n t o the r i v e r . A l a r g e low K zone at the toe of the slope may impede the flow of e f f l u e n t from the pond. The r e s u l t s of the s i m u l a t i o n s are presented i n Table 13. In a l l three r e f e r e n c e cases the b i a s i s the lowest f o r s i m u l a t i o n s based on sampling i n the region of the pond ( F ) . Each RMSE and the average RMSE are a l s o lowest f o r s i m u l a t i o n F. These r e s u l t s suggest that f o r the re f e r e n c e cases examined, key regions may e x i s t in the flow domain that are more important to sample than other r e g i o n s . For the re f e r e n c e cases with s m a l l to moderate standard d e v i a t i o n i n h y d r a u l i c c o n d u c t i v i t y and moderate to l a r g e i n t e g r a l s c a l e s , the region underneath the pond tends to be an important region to sample f o r reducing the u n c e r t a i n t y i n p r e d i c t i o n s of pond d i s c h a r g e . For some of the re f e r e n c e RESULTS / 120 100 -i Q l 3 80 -60 -40 20 / Q = 0.56 12 i r CO CO r-. CO d d Q 1 1 = 0 .93 ° 1 2 = 0 .82 ° 1 3 = 0.47 f = l F CNi d rr CO CO rr CVJ CO CO o <- co rf 00 lO cvj CO CO CO CO Is* i - CVJ cvj cvj Q (10 2 m 2 / s ) F i g . 23. Frequency d i s t r i b u t i o n of pond d i s c h a r g e from the a p r i o r i model \x=16m, (uy=-3, oy=0.5, X z=6m). Pond d i s c h a r g e from refere n c e cases 11 through 13 shown f o r comparison. RESULTS / 121 T A B L E 13. C o n d i t i o n a l S i m u l a t i o n s of Measurements i n Regions E through H f o r Reference Cases 11 through 13 RCN Y * cm av Qt Q Q-Qt s9 RMSE Region E 11 -3.26 0.50 -3.13 0.50 0.93 0.44 -0.49 0.24 0.54 0.55 12 -3.10 0.40 -3.13 6.39 0.82 0.38 -0.44 0.13 0.46 0.34 13 -3.17 0.48 -3.01 0.47' 0.47 0.54 0.08 0.24 0.25 0.44 Average 0.33** 0.20 0.42 Region F 11 -2.88 0.41 -2.86 0.40 0.93 0.81 -0.12 0.14 0.18 0.17 12 -2.81 0.38 -2.84 0.37 0.82 0.98 0.16 0.13 0.20 0.13 13 -3.32 0.38 -3.32 0.37 0.47 0.33 -0.13 0.04 0.14 0.13 Average 0.14 0.10 0.18 Region G 11 -2.93 0.49 -3.06 0.49 0.93 0.50 -0.43 0.25 0.49 0.50 12 -3.25 0.28 -3.23 0.28 0.82 0.29 -0.53 0.08 0.54 0.27 13 -2.60 0.46 -2.69 0.46 0.47 1.14 0.67 0.54 0.86 0.48 Average 0.54 0.29 0.63 Region H 11 -3.10 0.52 -3.01 0.51 0.93 0.G0 -0.33 0.32 0.46 0.53 12 -3.23 0.45 -3.15 0.44 0.82 0.38 -0.44 0.16 0.47 0.42 13 -3.84 0.43 -3.74 0.42 0.47 0.10 -0.37 0.04 0.37 0.41 Average 0.38 0.17 0.43 * Y, Sy for K in m/s; Q, RMSE in (m 2 / a X 10 - 2) ** E{\Q-Q,\] cases, however, the s p a t i a l v a r i a t i o n i n h y d r a u l i c c o n d u c t i v i t y may cause s i m u l a t i o n s based on samples from under the pond to le a d to l e s s c e r t a i n p r e d i c t i o n s of pond di s c h a r g e than other s i m u l a t i o n s . In those cases, RESULTS / 122 i d e n t i f y i n g the unknown p a t t e r n s of s p a t i a l v a r i a t i o n i n h y d r a u l i c c o n d u c t i v i t y may be more important than sampling beneath the pond. LOCATION FOR SECOND BOREHOLE Reference cases 11 through 13 R e s u l t s from s i m u l a t i o n s A through D suggest that an i n i t i a l borehole should be l o c a t e d below the cente r of the pond. S i m u l a t i o n s I through L are designed to i n v e s t i g a t e where a second borehole should be l o c a t e d given the l o c a t i o n of an i n i t i a l borehole (F igure 9). The d e c i s i o n on where to l o c a t e the second borehole i s based on the use of the data from both boreholes r a t h e r than on an a n a l y s i s of the data c o l l e c t e d i n the f i r s t b orehole. C o n d i t i o n a l s i m u l a t i o n s were run on ref e r e n c e cases with high standard d e v i a t i o n and small i n t e g r a l s c a l e s (11 through 13). The r e s u l t s are presented i n Table 14. In scheme I, both boreholes are l o c a t e d below the pond with the i n i t i a l borehole i n l o c a t i o n A below the center of the pond and the second borehole i n l o c a t i o n B below the edge of the pond. For a l l three r e f e r e n c e cases, the RMSE i s lower f o r scheme I than f o r scheme J , K, or L. In a d d i t i o n , the average RMSE i s lower f o r scheme I than f o r the other t h r e e . Caution must be used, however, i n i n t e r p r e t i n g the RESULTS / 123 TABLE 14. C o n d i t i o n a l S i m u l a t i o n s of Measurements i n Two Boreholes (I through L) f o r Reference Cases 11 through 13 RCN Qt Q Q-Qt Sv RMSE vq Scheme I 11 12 13 -2.75 0.50 -2.84 0.35 -3.11 0.46 -2.75 0.49 -2.82 0.35 -3.13 0.45 0.93 1.01 0.08 . 0.28 0.82 0.88 0.06 0.17 0.47 0.47 0.00 0.12 0.29 0.28 0.18 0.20 0.12 0.25 Average 0.05** 0.19 0.20 Scheme J 11 12 13 -2.41 0.58 -2.85 0.36 -3.12 0.32 -2.42 0.57 -2.83 0.36 -3.10 0.31 0.93 2.15 1.53 1.62 0.82 0.70 -0.11 0.21 0.47 0.37 -0.09 0.09 2.23 0.66 0.23 0.29 0.13 0.24 Average 0.58 0.63 0.85 Scheme K 11 12 13 -2.76 0.55 -3.13 0.43 -2.88 0.56 -2.80 0.53 -3.11 0.42 -2.90 0.56 0.93 1.02 0.09 0.57 0.82 0.43 -0.39 0.15 0.47 0.66 0.19 0.44 0.58 0.56 0.42 0.36 0.48 0.67 Average 0.23 0.39 0.49 Scheme L 11 12 13 -2.84 0.49 -3.12 0.59 -3.10 0.45 -2.84 0.48 -3.04 0.58 -3.10 0.44 0.93 0.88 -0.04 0.41 0.82 0.56 -0.26 0.33 0.47 0.39 -0.07 0.15 0.41 0.46 0.42 0.58 0.16 0.37 Average 0.13 0.30 0.33 * Y, Sy for A' in m/s; Q, RMSE in (m2/s x 10~2) **E\\Q-Qt\] average value of RMSE. Reference case 11 has an e x c e p t i o n a l l y high value of RMSE when scheme J i s used. T h i s is'because the second borehole i s l o c a t e d i n a hig h K region downstream of the pond and y i e l d s a l a r g e overestimate of RESULTS / 124 Y and an unusually l a r g e overestimate of Qt . The data f o r the i n d i v i d u a l r e f e r e n c e cases suggest that when the i n t e g r a l s c a l e s are sma l l , a second borehole l o c a t e d below the pond pro v i d e s more u s e f u l i n f o r m a t i o n f o r p r e d i c t i n g pond discharge than boreholes l o c a t e d elsewhere. SAMPLING IN STRATIFIED MEDIA Reference cases 14 through 17 are designed to i n v e s t i g a t e the e f f e c t i v e n e s s of d i f f e r e n t sampling s t r a t e g i e s i n s t r a t i f i e d media. These r e f e r e n c e cases are s t r o n g l y c o r r e l a t e d i n the h o r i z o n t a l d i r e c t i o n (XJC=160 m) and weakly c o r r e l a t e d i n the v e r t i c a l d i r e c t i o n (X z = 6 m). Due to the l a r g e standard d e v i a t i o n (a =0.5), a strong c o n t r a s t e x i s t s between the low K l a y e r s and the high K l a y e r s . Maps of l o g K, s t r e a m l i n e p l o t s , and f l u x p r o f i l e s f o r r e f e r e n c e cases 14 through 17 are shown in Figure 24. The contour i n t e r v a l f o r the four s t r e a m l i n e p l o t s i s 1.8 x 10" 3 m2/s . In re f e r e n c e case 14, the pond i s u n d e r l a i n by two 3 m t h i c k l a y e r s of low K with a 9 m t h i c k l a y e r of mid K in between. At depth, there i s an e x t e n s i v e high K l a y e r . A l l of the di s c h a r g e from the pond t r a v e l s downward i n t o the high K l a y e r and d i s c h a r g e s i n t o the r i v e r . Flux through the pond i s more uniform than i n the other r e f e r e n c e c a s e s . Qt i s the same as the average value from the a p r i o r i model (0.67 x 10" 2 m2/s, F i g u r e 25). RESULTS / 125 F i g . 24. Reference cases 14 through 17. (a) Map of l o g K. (b) Streamline p l o t with contour i n t e r v a l = 1.8 x 10' 3 m 2/s. (c) Flux p r o f i l e . ELEVATION, in EL£VATION. m ELEVATION, m ELEVATION, in X g o IS 8 0 8 8 0 8 8 921 RESULTS / 127 The pond in r e f e r e n c e case 15 i s d i r e c t l y u n d e r l a i n by a mid to high K l a y e r that o v e r l i e s a 6 to 12 m t h i c k l a y e r of low K. Some important d i f f e r e n c e s e x i s t between t h i s r e f e r e n c e case and r e f e r e n c e case 14. In r e f e r e n c e case 15, the low K l a y e r under the pond i s t h i c k e r and o v e r l i e s a high K l a y e r that i s not as e x t e n s i v e as the high K l a y e r i n r e f e r e n c e case 14. The high K l a y e r pinches out i n low K l a y e r s . F l u i d t r a v e l l i n g in the hi g h K l a y e r must a l s o pass through one of the low K l a y e r s . F i g u r e 26a pro v i d e s more d e t a i l of the s t r e a m l i n e p l o t , using a sm a l l e r contour i n t e r v a l . Streamlines are c l o s e together i n the high K l a y e r and d i v e r g e when the l a y e r pinches out. T h i s arrangement of h y d r a u l i c c o n d u c t i v i t y l a y e r s c o n s t r i c t s the flow and causes the lowest Qt of the three r e f e r e n c e cases. The l a y e r s below the pond i n r e f e r e n c e case 16 are, from top to bottom, a high K l a y e r , a mid to low K l a y e r , a high to mid K l a y e r , an e x t e n s i v e low K l a y e r , and a mid K l a y e r . Because the uppermost high K l a y e r i n t e r s e c t s the water t a b l e and i s u n d e r l a i n by a low K l a y e r , at l e a s t 40 percent of the e f f l u e n t from the pond e x i t s a c r o s s the water t a b l e j u s t downstream of the edge of the pond (Figure 26b). In the e x t e n s i v e low K l a y e r , the s t r e a m l i n e s d i v e r g e as e f f l u e n t t r i e s to bypass t h i s l a y e r and flow i n the mid K l a y e r s above and below the low K l a y e r . Qt f o r t h i s r e f e r e n c e case i s between that of r e f e r e n c e cases 14 and 15. RESULTS / 128 8 0 6 0 P. 4 0 2 0 0 ^ 1 5 Q / 14 i Q 0.67 X 10 2 m 2 / s — Sq = 0.33 X 10" 2 m 2 /s / °- Ql6 Q 1 4 0.62 X 10" 2 m 2 / s r ° 1 5 0.37 X 10"2 m 2 / s Q , 6 Q 1 7 — 0.67 0.71 X X 10" 2 m 2 /s 10" 2 m 2 / s - r — i — > I i i • i I — i 1 1— CD ' - C O ^ - C O T - C D O L O O C O O I O O I O O OTTfcOr^OiOOJCOiO'<- CM d 0 0 0 0 0 * - i - T - T - T - T - i T - ' c N i c J c v i Q ( 1 0 _ 2 m 2 / s ) F i g . 25. Frequency d i s t r i b u t i o n of pond d i s c h a r g e from the a p r i o r i model (y =-3, a^=0.5, Xx=160m, Xz=6m). Pond d i s c h a r g e from r e f e r e n c e cases 14 through 17 shown f o r comparison. RESULTS / 129 a . o ^ 1 1 1 1 1 1 1 1 1 r b . 0 H 1 1 1 1 1 1 1 1 i i 0 30 60 90 120 150 180 210 240 270 300 HORIZONTAL DISTANCE, metres F i g . 26. (a) Streamline p l o t f o r r e f e r e n c e case 15 with contour i n t e r v a l 4.75 x 10'" m 2/s. (b) Streamline p l o t f o r r e f e r e n c e case 16 with contour i n t e r v a l 8.0 x 10"" m 2/s. RESULTS / 130 The highest Qt i s from r e f e r e n c e case 17. There are b a s i c a l l y three l a y e r s ; two high K l a y e r s separated by a low K l a y e r . As in re f e r e n c e case 14, a l l of the e f f l u e n t from the pond e v e n t u a l l y t r a v e l s down i n t o the b a s a l high K l a y e r and d i s c h a r g e s i n t o the r i v e r . Table 15 l i s t s the r e s u l t s of a p p l y i n g sampling schemes M through O to the s t r a t i f i e d r e f e r e n c e cases. Scheme M has three measurements i n each of s i x boreholes spaced evenly throughout the flow domain (Figure 10). The purpose of scheme N and 0 i s to determine whether more measurements should be taken than 18 (scheme M) and whether the number of boreholes should be i n c r e a s e d (scheme N) or whether the number of samples i n each borehole should be in c r e a s e d (scheme 0 ) . In r e f e r e n c e case 14, the s i m u l a t i o n s based on the 18 measurements from sampling scheme M y i e l d e d a lower RMSE than e i t h e r scheme N or 0 with 30 measurements. R e c a l l that sampling schemes N and 0 i n c l u d e the 18 measurements of scheme M. Adding 4 a d d i t i o n a l boreholes (scheme N) f a i l e d to change the RMSE by much, p a r t l y because of the l a r g e h o r i z o n t a l i n t e g r a l s c a l e . Adding 2 more measurements to each of the s i x boreholes of scheme M l e d to a l e s s c e r t a i n p r e d i c t i o n of pond d i s c h a r g e . N e i t h e r sampling schemes M or N c o l l e c t data from e i t h e r of the two low K l a y e r s underneath the pond, which leads to overestimates of pond RESULTS / 131 TABLE 15. C o n d i t i o n a l S i m u l a t i o n s of Measurements using Schemes M through 0 f o r Reference Cases 14 through 17 RCN Y * cm °y Yc Qt Q Q-Qt Sq RMSE Scheme M 14 -2.70 0.43 -2.74 0.42 0.62 1.06 0.44 0.50 0.67 0.47 15 -3.31 0.35 -3.3S 0.33 0.37 0.28 -0.09 0.11 0.14 0.40 16 -3.25 0.28 -3.27 0.25 0.67 0.27 -0.40 0.07 0.40 0.24 17 -2.86 0.48 -2.80 0.49 0.71 1.28 0.57 0.62 0.84 0.49 Average 0.37** 0.33 0.51 0.40 Scheme N 14 -2.72 0.37 -2.67 0.36 0.62 1.14 0.53 0.43 0.68 0.37 15 -3.23 0.35 -3.31 0.32 0.37 0.29 -0.08 0.10 0.12 0.34 16 -3.28 0.28 -3.28 0.27 0.67 0.29 -0.38 0.08 0.39 0.29 -2.89 0.43 -2.80 0.47 0.71 1.24 0.52 0.58 0.78 0.47 Average 0.38 0.30 0.49 0.37 Scheme 0 14 -2.69 0.50 -2.77 0.49 0.62 1.06 0.44 0.62 0.76 0.59 15 -3.29 0.38 -3.33 0.35 0.37 0.29 -0.08 0.11 0.13 0.37 16 -3.32 0.33 -3.32 0.30 0.67 0.26 -0.41 0.09 0.42 0.34 17 -2.98 0.46 -2.81 0.48 0.71 1.20 0.48 0.57 0.75 0.47 Average 0.35 0.35 0.52 0.44 * Y, Su for A' in m/s; Q, RMSE in (m2/s x H T 2 ) **E[\Q-Q,\] d i s c h a r g e . Although scheme 0 c o n t a i n s measurements from these low K l a y e r s i t a l s o c o n t a i n s a d d i t i o n a l measurements, from the exte n s i v e high K l a y e r , which leads to an overestimate of Qt as w e l l . In r e f e r e n c e case 15, the RMSE i s small f o r the three s i m u l a t i o n s . A l l of the sampling schemes miss the t h i n high RESULTS / 132 K l a y e r and l e a d to s l i g h t underestimates of Q{. Because the 18 measurements i n scheme M l e a d to a low RMSE, i n c r e a s i n g the number of measurements (scheme N and 0) does l i t t l e to improve an alr e a d y low RMSE. A l l three sampling schemes l e a d to lower RMSEs f o r r e f e r e n c e case 15 than f o r the other three r e f e r e n c e cases. Reference case 15 has a smoother f l u x p r o f i l e than the other s t r a t i f i e d r e f e r e n c e cases. It may be e a s i e r to p r e d i c t pond discharge when f l u x a c r o s s the water t a b l e i s small and v a r i e s only s l i g h t l y . In r e f e r e n c e case 16, l i t t l e change i s observed i n the RMSE between schemes. The p r e d i c t i o n s of pond discharge are l e s s c e r t a i n f o r r e f e r e n c e case 16 than f o r 15. None of the three sampling schemes c o l l e c t samples i n the hi g h K l a y e r d i r e c t l y under the pond and a l l three l e a d to underestimates of Qt . Both schemes M and N miss the e x t e n s i v e low K l a y e r at depth. Scheme 0 c o n t a i n s measurements i n t h i s low K l a y e r , which leads to a lower estimate of Y , an underestimate of Q. of a higher magnitude, a higher S , and i q a s l i g h t l y higher RMSE. The hi g h e s t RMSEs f o r sampling schemes N through 0 are from r e f e r e n c e case 17. A l l three sampling schemes l e a d to overestimates of Qt , with scheme M y i e l d i n g the hi g h e s t b i a s and RMSE, and scheme 0 y i e l d i n g the lowest. Compared to the other sampling schemes, the percentage of measurements i n the h i g h K l a y e r s i s g r e a t e s t and the percentage of RESULTS / 133 measurements i n the low K l a y e r s i s sm a l l e s t f o r scheme M. The op p o s i t e i s true f o r scheme 0 with the lowest RMSE. In two of the four r e f e r e n c e cases, the RMSE changes only s l i g h t l y f o r the three sampling schemes. In re f e r e n c e case 14, i t i s more e f f e c t i v e to take 30 measurements i n ten boreholes (scheme N) than i n s i x boreholes (scheme 0 ) . The cos t of scheme N, however, i s much hig h e r . In re f e r e n c e case 17, scheme 0 y i e l d s a lower RMSE than scheme N, because the v e r t i c a l l y spaced data i n scheme O d e l i n e a t e s the l a y e r s b e t t e r . Because the h o r i z o n t a l i n t e g r a l s c a l e i s l a r g e , the d i f f e r e n c e i n RMSE between s i m u l a t i o n s using s i x boreholes with three measurements i n each and s i m u l a t i o n s using ten boreholes with three measurements i n each i s g e n e r a l l y only s l i g h t . The co s t of the ten boreholes (N) i s much gr e a t e r than the co s t of the s i x boreholes (M). The three sampling schemes d e s c r i b e d i n t h i s s e c t i o n appear most e f f e c t i v e i n reducing p r e d i c t i o n u n c e r t a i n t y i n r e f e r e n c e case 15 and l e a s t e f f e c t i v e i n r e f e r e n c e case 17. In summary, i n c r e a s i n g the number of measurements to 30 from 18 by e i t h e r d r i l l i n g new boreholes or t a k i n g more measurements i n each borehole w i l l not n e c e s s a r i l y l e a d to s i g n i f i c a n t r e d u c t i o n s i n p r e d i c t i o n u n c e r t a i n t y i n s t r a t i f i e d media. T h i s may be p a r t l y due to the strong h o r i z o n t a l c o r r e l a t i o n . Although the u n c e r t a i n t y i n K i s reduced over a l a r g e d i s t a n c e i n s t r a t i f i e d media, RESULTS / 134 c o n s i d e r a b l e u n c e r t a i n t y i n d i s c h a r g e p r e d i c t i o n s can s t i l l e x i s t even with a r e l a t i v e l y l a r g e number of measurements. When the h y d r a u l i c c o n d u c t i v i t y measurements are spaced evenly through the flow domain, important shallow l a y e r s may be missed that have a l a r g e i n f l u e n c e on pond d i s c h a r g e . MULTIPLE BOREHOLES The l a s t set of s i m u l a t i o n s i s designed to i n v e s t i g a t e the e f f e c t i v e n e s s of m u l t i p l e boreholes i n reducing p r e d i c t i o n u n c e r t a i n t y . I f m u l t i p l e boreholes are to be d r i l l e d , at what p o i n t does the i n c r e a s e i n c o s t of a d d i t i o n a l boreholes outweigh the gain i n information? To address t h i s q u e s t i o n l e t us now r e t u r n to r e f e r e n c e cases 6 through 10 using sampling schemes P through R. Scheme P has four boreholes, scheme Q has seven boreholes and scheme R has 13 boreholes. Because the p r e v i o u s sampling schemes (M, N, and 0) f r e q u e n t l y missed important l a y e r s d i r e c t l y underneath the pond, the measurements i n schemes P through R are p l a c e d c l o s e r to the pond and water t a b l e than to the b a s a l boundary. Table 16 l i s t s the r e s u l t s of the c o n d i t i o n a l s i m u l a t i o n s . In r e f e r e n c e case 6, the RMSE i s very low f o r scheme P with the s m a l l e s t number of boreholes. No advantage i s gained by d r i l l i n g a d d i t i o n a l b o r eholes. In r e f e r e n c e case 7, the RMSE i s lowest f o r 4 boreholes and much higher f o r 7 RESULTS / 135 TABLE 16. C o n d i t i o n a l S i m u l a t i o n s of Measurements i n M u l t i p l e Boreholes P through R f o r Reference Cases 6 through 10 RCN Y * 1 m cm Ye CC Qt Q Q-Qt Sq RMSE vq Scheme P 6 -3.02 0.31 -3.04 0.3Q 0.3S 0.39 -0.01 0.13 0.13 0.32 7 -2.72 0.31 -2.76 0.29 1.45 LOG -0.38 0.35 0.52 0.33 8 -2.99 0.40 -2.95 0.38 0.55 0.69 0.15 0.29 0.33 0.42 9 -2.87 0.44 -2.84 0.41 0.41 0.78 0.37 0.33 0.50 0.42 10 -3.17 0.39 -3.14 0.38 0.42 0.41 -0.01 0.17 0.17 0.41 Average 0.33 0.38 Standard Deviation 0.18 0.05 Scheme Q 6 -3.06 0.43 -3.07 0.42 0.38 0.29 -0.09 0.11 0.14 0.38 7 -2.66 0.42 -2.G9 0.41 1.45 1.72 0.27 0.67 0.72 0.39 8 -2.90 0.40 -2.89 0.38 0.55 0.67 0.12 0.22 0.25 0.34 9 -3.01 0.42 -2.99 0.39 0.41 0.39 -0.03 0.13 0.13 0.32 10 -3.15 0.40 -3.11 0.39 0.42 0.53 0.10 0.17 0.20 0.33 Average 0.29 0.35 Standard Deviation 0.25 0.03 Scheme R. 6 -3.08 0.41 -3.08 0.41 0.38 0.29 -0.09 0.10 0.14 0.35 7 -2.61 0.47 -2.62 0.46 1.45 1.72 0.27 0.65 0.71 0.38 8 -2.98 0.43 -2.98 0.41 0.55 0.57 0.02 0.17 0.17 0.31 9 -3.06 0.44 -3.02 0.42 0.41 0.40 -0.02 0.13 0.13 0.33 10 -3.11 0.38 -3.13 0.35 0.42 0.49 0.06 0.15 0.17 0.32 Average 0.26 0.34 Standard Deviation 0.25 0.03 * Y, Sy for A' in m/s; Q, RMSE in (m2/s x 10~2) **E\\Q-Qt\] and 13 boreholes. Even though the magnitude of the b i a s i s l a r g e r f o r the smaller number of boreholes, the standard RESULTS / 136 d e v i a t i o n of the measurements i s much s m a l l e r , l e a d i n g to the lower RMSE. In r e f e r e n c e case 8, each a d d i t i o n a l set of boreholes reduces the u n c e r t a i n t y i n p r e d i c t i o n s of Qt , by reducing both the b i a s and S^. In r e f e r e n c e case 9, the RMSE i s g r e a t l y reduced i n going from 4 to 7 boreholes, but remains the same f o r 7 and 13 boreholes. T h i s may be due to the h o r i z o n t a l c o r r e l a t i o n i n h y d r a u l i c c o n d u c t i v i t y . When the spacing of the boreholes approaches the h o r i z o n t a l i n t e g r a l s c a l e , i n c r e a s i n g the number of boreholes, may f a i l to p r o v i d e new i n f o r m a t i o n . In r e f e r e n c e case 10, the RMSE i s low f o r 4 boreholes and does not improve by adding a d d i t i o n a l boreholes. To summarize the r e s u l t s , the p r e d i c t i o n u n c e r t a i n t y i n e s t i m a t i n g pond di s c h a r g e decreases as the number of boreholes i n c r e a s e s f o r r e f e r e n c e cases 8 and 9. For re f e r e n c e cases 6, 7, and 10, there i s no advantage to be gained by c o l l e c t i n g data i n more than 4 boreholes. T h i s may be due i n part to the r e l a t i v e l y l a r g e h o r i z o n t a l i n t e g r a l s c a l e which reduces the u n c e r t a i n t y i n h y d r a u l i c c o n d u c t i v i t y over a l a r g e d i s t a n c e . The average RMSE from the f i v e r e f e r e n c e cases decreases from 4 boreholes to 7 boreholes and f u r t h e r decreases from 7 to 13 boreholes. Because of the small sample s i z e , however, these r e s u l t s must be viewed with c a u t i o n . RESULTS / 137 OBJECTIVE FUNCTION The r e s u l t s presented thus f a r compare sampling s t r a t e g i e s on the b a s i s of the RMSE, a measure of the p r e d i c t i o n u n c e r t a i n t y . We now wish to i n c o r p o r a t e economic c o n s i d e r a t i o n s i n t o comparisons of the sampling s t r a t e g i e s . An o b j e c t i v e f u n c t i o n (14) c o n s i d e r s both the c o s t of c o l l e c t i n g data, C^, and the l o s s e s a s s o c i a t e d with poor p r e d i c t i o n s , C^. R e c a l l that the l o s s f u n c t i o n i s the product of a co s t c o e f f i c i e n t , and the RMSE. i s the product of the width of the pond ( L ^ ) , the time h o r i z o n (Tn), and the c o s t of a l t e r n a t i v e treatment, Ct . R e a l i s t i c c o s t s f o r a l t e r n a t i v e treatment were c a l c u l a t e d f o r the co s t of c o l o r removal from pulp m i l l e f f l u e n t u s i n g the Skookumchuck pulp m i l l i n Cranbrook, B r i t i s h Columbia [Swaney and Jackson, 1983], The m i l l produces 450 a i r dry tons of bleached k r a f t (a.d.t.bk) d a i l y from 14.4 m i l l i o n U. S. g a l l o n s of e f f l u e n t . The c o s t of c o l o r removal i s estimated at $ 1.08/a.d.t.bk f o r r a p i d i n f i l t r a t i o n and ranges from $2.33 to $ 12.05/a.d.t.bk f o r c o n v e n t i o n a l methods. These f i g u r e s convert t o $8.92 x l0" 3/m 3 f o r treatment u s i n g r a p i d i n f i l t r a t i o n ponds, and $1.92 x l0- 2/m 3 to $9.95 x l0" 2/m 3 f o r c o n v e n t i o n a l treatment. Using a pond width of 10 m and a time h o r i z o n of 5 yea r s , C ; i s $1.41 x I0 7s/m 2 f o r r a p i d i n f i l t r a t i o n and from $3.03 x I0 7s/m 2 to $1.57 x I0 8s/m 2 f o r RESULTS / 138 c o n v e n t i o n a l methods. Table 17 presents the c o s t of sampling (C^), RMSE, value of the l o s s f u n c t i o n (C^), and value of the o b j e c t i v e f u n c t i o n (Z) f o r each of the c o n d i t i o n a l s i m u l a t i o n s . Because of the h y p o t h e t i c a l nature of t h i s study, an a r b i t r a r y and minimal value of ($1.30 x I0 7s/m 2) i s used. The l a s t column of Table 17 l i s t s the s e n s i t i v i t y of the r e s u l t s to changes in the value of C^. The cost of c o l l e c t i n g data i n the i n i t i a l boreholes A, B, C, and D i s $8,315, $8,315, $10,714 and $11,996. The cost of measurements i n boreholes A and B i s l e s s than i n C or D because the former have a higher percentage of unsaturated samples. No attempt i s made to i n c o r p o r a t e measurement e r r o r or the d i f f e r e n c e i n r e l i a b i l i t y between measuring K from unsaturated samples i n the l a b o r a t o r y and measuring s a t u r a t e d K i n - s i t u . An a s t e r i s k (*) i n the l a s t column of Table 17 i n d i c a t e s that the corresponding sampling s t r a t e g y i s always the best s t r a t e g y of the group of s t r a t e g i e s t e s t e d , r e g a r d l e s s of the value of . T h i s occurs when both the RMSE and the c o s t of sampling, C$, are lower f o r one s t r a t e g y than f o r the other s t r a t e g i e s . The minimum value of the o b j e c t i v e f u n c t i o n i s then independent of the value of C[. A blank entry i n t h i s column i n d i c a t e s that the c o r r e s p o n d i n g sampling s t r a t e g y i s never the best s t r a t e g y , RESULTS / 139 TABLE 17. O b j e c t i v e F unction f o r the C o n d i t i o n a l Simulations Reference case 1 Scheme Cs RMSE CL Z Best design A 8,300 0.17 22,100 30,400 * B 8,300 C.25 32,500 40,800 C 10,700 0.18 23,400 34,100 D 12,000 , 0.18 23,400 35,400 E 34,200 0.42 54,600 88,800 . F 23,800 0.26 33,800 57,600 Gt < 2.o7 x 10° G 26,200 0.16 20,800 47,000 2.37 x 106 < Ci < 2.43 x 107 H 28,600 0.15 19,500 48,100 C ( > 2.43 x 107 Reference case 2 A 8,300 0.04 5,200 13,500 * B 8,300 0.16 20,800 29,100 C 10,700 0.08 10,400 21,100 D 12,000 0.28 36,400 48,400 E 34,200 0.11 14,300 48,500 F 23,800 0.06 7,800 31,600 * G 26,200 0.06 7,800 34,000 H 28,600 0.25 32,500 . 61,100 Reference case 3 A 8,300 0.07 9,100 17,400 * B 8,300 0.09 11,700 20,000 C 10,700 0.26 33,800 44,500 D 12,000 0.10 13,000 25,000 -E 34,200 0.18 23,400 57,600 F 23,800 0.07 9,100 32,900 * G 26.200 0.25 32,500 58,700 H 28,600 0.08 10,400 39,000 Reference case 4 A 8,300 0.06 7,800 16,100 * B 8,300 0.11 14,300 22,600 C 10,700 0.07 9,100 19,800 D 12,000 0.18 23,400 35,400 E 34,200 0.14 18,200 52,400 F 23,800 0.07 9,100 32,900 * G 26,200 0.07 9,100 35,300 H 28,600 0.31 40,300 68,000 Cs, CL, Z in $; RMSE in (m2/s * 10~2). Based on Ct = $8.26 x 10" 5/™ 3, Lw = 10m, T = 5 years; C, = $1.30 x 107«/m2. RESULTS / 140 TABLE 17. (continued) Reference case 5 Scheme Gs RMSE CL Z Best design A 8,300 0.16 20,800 29,100 B 8,300 0.06 7,800 16,100 C 10,700 0.25 32,500 43,200 D 12,000 0.19 • 24,700 36,700 E 34,200 0.22 28,600 62,800 F 23,800 0.11 14,300 38,100 Ct < 2.40 x 107 G 26,200 0.23 29,900 56,100 H 28,600 0.09 11,700 40,300 Ci > 2.40 x 107 Reference case 6 E 34,200 0.16 20,800 55,000 Ci > 3.45 x 107 F 23,800 0.19 24,700 48,500 Ci < 3.45 x 107 G 26,200 0.70 91,000 117,200 H 28,600 0.65 84,500 113,100 P 20,700 0.13 16,900 37,600 * Q 34,200 0.14 18,200 52,400 R 61,600 0.14 18,200 79,800 Reference case 7 E 34,200 1.26 163,800 198,000 F 23,800 0.40 52,000 75,800 * G 26,200 0.66 85,800 112,800 H 28,600 0.86 111,800 140,400 P 20,700 0.52 67,600 88,300 * Q 34,200 0.72 93,600 127,800 R 61,600 0.71 92,300 153,900 Reference case 8 E 34,200 0.27 35,100 69,300 F 23,800 0.29 37,700 61,500 Ci < 4.80 x 106 G 26,200 0.32 41,600 67,800 H 28;G00 0.19 24,700 53,300 Ci > 4.80 x 107 P 20,700 0.33 42,900 63,600 Ci < 1.70 x 107 Q 34,200 0.25 32,500 66,700 1.70 x 107 < Ci < 3.42 x 107 R 61,600 0.17 22,100 83,700 Ci > 3.42 x 107 Cs, CL, Z in $; RMSE in ( m 2 / * x 10~2). Based on Ct = $8.26 x 1 0 _ 5 / m s , LW = 10m, T = 5 years; Ct = $1.30 X 1 0 7 « / m 2 . RESULTS / 141 TABLE 17. (continued) Reference case 9 Scheme CS RMSE CL Z i Best design E 34,200 0.26 33,800 68,000 F 23,800 0.15 19,500 43,300 * G 26,200 0.25 32,500 58,700 H 28,600 0.25 32^ 500 61,100 P 20,700 0.50 65,000 85,700 Ci < 3.66 x 106 Q 34,200 0.13 16,900 51,100 Ci > 3.66 x 106 R 61,600 0.13 16,900 78,500 Reference case 10 E 34,200 0.21 27,300 61,500 F 23,800 0.11 14,300 38,100 * G 26,200 0.33 42,900 69,100 H 28,000 0.22 28.600 57,200 P 20,700 0.17 22,100 42,800 * Q 34,200 0.20 26,000 60,200 R 61,600 0.17 22,100 83,700 Reference case 11 E 34,200 0.54 70,200 104,400 F 23,800 0.18 23,400 47,200 * G 26,200 0.49 63,700 89,900 H 28,000 0.46 59,800 88,400 I 11,000 0.29 37,700 48,700 * J 11,800 2.23 289,900 301,700 K 12,600 0.58 75,400 88,000 L 13,400 0.41 53,300 66,700 Reference case 12 E 34,200 0.46 59,800 94,000 F 23,800 0.20 26,000 49,800 * G 20,200 0.54 70,200 96,400 H 28,600 0.47 61,100 89,700 I 11,000 0.18 23,400 34,400 * J 11,800 0.23 29,900 41,700 K 12,600 0.42 54,600 67,200 L 13,400 0.42 54,000 68,000 C s , Ct, Z in $; RMSE in (m2/* x 10 - 2). Based on C, = $8.26 x 10-s/ms, Lw = 10m, T = 5 years; C ( = $1.30 x 107«/m2. RESULTS / 142 TABLE 17. (continued) Reference case 13 Scheme Cs RMSE CL Z Best design E 34,200 0.25 32,500 66,700 F 23,800 0.14 18,200 42,000 * G 26,200 0.86 111,800 138,000 H 28,600 0.37 48,100 76,700 I 11,000 0.12 15,600 26,600 * J 11,800 0.13 16,900 28,700 K 12,600 0.48 62,400 75,000 L 13,400 0.16 , 20,800 34,200 Reference case 14 M 33,200 0.67 87,100 120,300 * N 53,100 0.68 88,400 141,500 O 38,300 0.76 98,800 137,100 Reference case 15 M 33,200 0.14 18,200 51,400 Ci < 5.09 X 107 N 53,100 0.12 15,600 68,700 Ci > 1.48 X 108 O 38,300 0.13 16,900 55,200 5.09 x 107 < Ci < 1.48 x 108 Reference case 16 M 33,200 0.40 52,000 85,200 Ci < 1.99 x 108 N 53,100 0.39 50,700 103,800 Ci > 1.99 x 108 O 38,300 0.42 54,600 92,900 Reference case 17 M 33,200 0.84 1.09 x 107 1.10 x 107 Ci < 5.66 X iO 6 N 53,100 0.78 1.01 x 107 1.02 x 107 O 38,300 0.75 9.75 x 106 9.79 x 106 Ci > 5.66 x 106 Cs, CL, Z in $; RMSE in (m2/g x 10~2). Based on Ct = $8.26 x lO"*/™8, Lw = 10m, T = 5 years; (7, = $1.30 x 107a/m2. RESULTS / 143 r e g a r d l e s s of the value of C[. Consider r e f e r e n c e case 1. The minimum value of the o b j e c t i v e f u n c t i o n i s from borehole A l o c a t e d below the middle of the pond when i s used as input. In r e f e r e n c e case 1, borehole A i s always the best sampling s t r a t e g y of A through D because i t has the lowest RMSE and the sampling c o s t s of the other s t r a t e g i e s are the same or h i g h e r . For r e f e r e n c e cases 2, 3, and 4, borehole A, l o c a t e d below the middle of the pond, always y i e l d s the minimum value of Z f o r the s t r a t e g i e s t e s t e d , r e g a r d l e s s of the value of C[. For r e f e r e n c e case 5, borehole B i s the best s t r a t e g y of the fou r . In summary, sampling below the pond tends to be the best s t r a t e g y f o r the cases t e s t e d . Borehole A, l o c a t e d below the middle of the pond, i s p r e f e r r e d f o r four of the ref e r e n c e cases and borehole B, below the edge of the pond, i s p r e f e r r e d f o r one of the ref e r e n c e cases. For sampling i n regions E through H f o r re f e r e n c e case 1, the best s t r a t e g y depends upon the co s t c o e f f i c i e n t , . Sampling i n region F (below the pond) i s the best s t r a t e g y when i s l e s s than $2.37 x I0 6s/m 2, sampling i n r e g i o n H (near the r i g h t hand s i d e ) i s p r e f e r r e d when i s g r e a t e r than $2.43 x I0 7s/m 2, and sampling i n region G (mid domain) i s p r e f e r r e d when Cj i s between these two v a l u e s . Sampling in r e g i o n E (lower l e f t ) i s never p r e f e r r e d i n t h i s case because i t y i e l d s the hi g h e s t RMSE and has the h i g h e s t RESULTS / 144 sampling c o s t s . Region F i s the best p l a c e to sample f o r re f e r e n c e cases 2, 3, and 4 and i s p r e f e r r e d in r e f e r e n c e case 5 f o r small Cj . For ay=0.5, region F i s p r e f e r r e d f o r re f e r e n c e cases 7, 9, and 10 and f o r re f e r e n c e cases 6 and 8, when Cf i s s m a l l . When the degree of s p a t i a l c o r r e l a t i o n i s s m a l l , F i s p r e f e r r e d f o r the three r e f e r e n c e cases t e s t e d (11 through 13). These r e s u l t s i n d i c a t e that f o r r e f e r e n c e cases 1 through 13, sampling in the region of the pond i s g e n e r a l l y p r e f e r r e d f o r p r e d i c t i n g pond d i s c h a r g e . In a d d r e s s i n g the qu e s t i o n of how many boreholes to d r i l l when the s p a t i a l c o r r e l a t i o n i s l a r g e , but not l a r g e enough to c r e a t e e x t e n s i v e l a y e r s , 4 boreholes (scheme P) i s best f o r ref e r e n c e cases 6, 7, and 10. For r e f e r e n c e cases 8 and 9, 4 boreholes i s b e t t e r than 7 or 13 when i s s m a l l . These r e s u l t s suggest that i f boreholes are spaced evenly throughout the flow domain, four boreholes are g e n e r a l l y s u f f i c i e n t f o r the r e f e r e n c e cases c o n s i d e r e d . I n c r e a s i n g the number of boreholes beyond four i s not worth the added expense i n most of the cases c o n s i d e r e d u n l e s s the cost c o e f f i c i e n t , C[, i s h i g h . When the h o r i z o n t a l i n t e g r a l s c a l e i s l a r g e enough to c r e a t e e x t e n s i v e l a y e r s , c o l l e c t i n g 30 samples from ten boreholes ( s t r a t e g y N) y i e l d s more c e r t a i n p r e d i c t i o n s of pond d i s c h a r g e i n three of the four r e f e r e n c e cases than RESULTS / 145 c o l l e c t i n g 30 samples from s i x boreholes ( s t r a t e g y 0 ) . For the value of used ($1.30 x I0 7s/m 2), however, the added co s t of the ex t r a boreholes i s not o f f s e t by an e q u i v a l e n t decrease i n the value of the l o s s f u n c t i o n , r e s u l t i n g i n values of Z that are smaller f o r scheme 0 than scheme N. Note that these r e s u l t s depend s t r o n g l y upon the value of Cj . For h i g h values of s t r a t e g y N would be p r e f e r r e d over s t r a t e g y 0 f o r three of the r e f e r e n c e cases, because of the small e r values of RMSE. C o l l e c t i n g 18 samples i n s i x boreholes (M) i s p r e f e r r e d f o r r e f e r e n c e case 14 and f o r small v a l u e s of f o r the three other r e f e r e n c e cases because of the lower sampling c o s t . More s i m u l a t i o n s would be necessary to determine whether s t r a t e g i e s M, N, or 0 tend to be p r e f e r r e d . When data are c o l l e c t e d from two boreholes and the s p a t i a l c o r r e l a t i o n i s s m a l l , l o c a t i n g both boreholes below the pond (scheme I) i s p r e f e r r e d f o r re f e r e n c e cases 11, 12, and 13. In a l l three cases, the r e s u l t s are independent of the value of . Strat e g y I leads to lower v a l u e s of RMSE and i s l e s s expensive than s t r a t e g i e s J , K, and L. It i s important to emphasize that the r e s u l t s from these s i m u l a t i o n s apply to s p e c i f i c r e f e r e n c e cases. Using the framework introduced here, a l a r g e number of s i m u l a t i o n s are necessary to draw c o n c l u s i o n s on optimal sampling s t r a t e g i e s . LIMITATIONS In i n t e r p r e t i n g the r e s u l t s of t h i s study, i t i s important to keep in mind the l i m i t a t i o n s . T h i s study i s conceptual i n nature and can not be a p p l i e d d i r e c t l y to a f i e l d study. The framework developed here can be used to i n v e s t i g a t e the e f f e c t i v e n e s s of sampling s t r a t e g i e s as they apply to s p e c i f i c ensemble parameters for a p a r t i c u l a r boundary value problem. The o b j e c t i v e f u n c t i o n depends upon knowing the t r u e value of d i s c h a r g e from the pond, Qt . In a f i e l d study, Qt would not be known durin g the pre-pond sampling phase. The framework i s u s e f u l as a f i r s t step i n e s t a b l i s h i n g g u i d e l i n e s f o r the design of sampling programs to p r e d i c t groundwater d i s c h a r g e . Another l i m i t a t i o n of t h i s study i s the use of a f i x e d water t a b l e f o r each s i m u l a t i o n . The f l u x p r o f i l e s and s t r e a m l i n e p l o t s from the r e f e r e n c e cases suggest that the p o s i t i o n of the water t a b l e would depend upon the s p a t i a l arrangement of h y d r a u l i c c o n d u c t i v i t y i n each case. By f i x i n g the water t a b l e , f l u i d i s allowed to e x i t the flow domain to a v o i d low h y d r a u l i c c o n d u c t i v i t y v a l u e s near the water t a b l e . T h i s may a f f e c t the c a l c u l a t e d value of d i s c h a r g e from the pond. In r e a l i t y , the water t a b l e would r i s e or f a l l i n response to the p a t t e r n s of s p a t i a l v a r i a b i l i t y i n K, and the e f f e c t of these low K zones near the water t a b l e would l i k e l y be evident in the value of 146 LIMITATIONS / 147 d i s c h a r g e from the pond. The use of a f r e e - s u r f a c e model would allow the water t a b l e to a d j u s t v e r t i c a l l y . Future s t u d i e s i n v e s t i g a t i n g s p a t i a l v a r i a b i l i t y on h i l l s l o p e s should combine s t o c h a s t i c methods with a f r e e s u r f a c e model. As shown in t h i s study, the c h o i c e of a sampling scheme i s dependent upon the s p a t i a l arrangement of h y d r a u l i c c o n d u c t i v i t y . To e s t a b l i s h g u i d e l i n e s f o r the design of sampling programs, a very l a r g e number of s i m u l a t i o n s would be r e q u i r e d . T h i s would be a monumental task, beyond the scope of t h i s study. Even given a l a r g e number of s i m u l a t i o n s with which to c a l c u l a t e the p r o b a b i l i t y d i s t r i b u t i o n of the RMSE f o r d i f f e r e n t s t r a t e g i e s , the r e s u l t s would only apply to the p a r t i c u l a r ensemble parameters and the p a r t i c u l a r boundary value problem chosen. T h e r e f o r e , i n i n t e r p r e t i n g the r e s u l t s from t h i s study i t i s e s s e n t i a l to bear in mind that they are s p e c i f i c to a small number of refe r e n c e cases. Another l i m i t a t i o n i s that the h o r i z o n t a l and v e r t i c a l i n t e g r a l s c a l e s are assumed known. A l a r g e number of K measurements are necessary to p r o p e r l y estimate these parameters. The u n c e r t a i n t y i n e s t i m a t i n g the i n t e g r a l s c a l e s are not i n c o r p o r a t e d i n t o the a n a l y s i s . F i n a l l y , t h i s study addresses the problem of network design i n only two dimensions. For an a c t u a l sampling program, the l o c a t i o n of K measurements in three dimensions LIMITATIONS / 148 would have to be c o n s i d e r e d . In a d d i t i o n , the p a t t e r n s of groundwater flow and di s c h a r g e would be expected to vary i n response to h e t e r o g e n e i t i e s in the t h i r d dimension. SUMMARY CONCLUSIONS 1. A framework was developed to a s s i g n economic worth to measurements of h y d r a u l i c c o n d u c t i v i t y . Using t h i s framework, the t r a d e o f f between the c o s t of data c o l l e c t i o n and the value of in f o r m a t i o n gained can be examined. 2. Network design f o r c o l l e c t i n g data that w i l l be used as input to a p r e d i c t i v e model must c o n s i d e r the r e d u c t i o n i n output u n c e r t a i n t y . Reducing the v a r i a n c e of e s t i m a t i o n f o r input parameters does not n e c e s s a r i l y l e a d to a r e d u c t i o n i n p r e d i c t i o n u n c e r t a i n t y . 3. For the problem c o n s i d e r e d , no optimal network design c o u l d be determined because a l a r g e number of r e f e r e n c e cases are necessary. P r e d i c t i o n u n c e r t a i n t y was s e n s i t i v e to both the s t r u c t u r e of the h e t e r o g e n e i t i e s and the l o c a t i o n of measurements. 4. The average estimate of di s c h a r g e i s i n s e n s i t i v e to the input standard d e v i a t i o n . The average estimate of pond di s c h a r g e i s more s e n s i t i v e to the average value of the measurements, than to the l o c a t i o n of the measurements. Measurements that have an average value c l o s e to that of the r e f e r e n c e case, however, may l e a d to poor p r e d i c t i o n s of di s c h a r g e because the s p e c i f i c 149 SUMMARY / 150 arrangement of h y d r a u l i c c o n d u c t i v i t y values can r e s u l t in v a l u e s of Qt that are lower or higher than the average value from the a p r i o r i model. T h e r e f o r e , sampling schemes which l e a d to good estimates of the ensemble mean and standard d e v i a t i o n w i l l not n e c e s s a r i l y lead to good p r e d i c t i o n s of d i s c h a r g e . The goal of sampling s t r a t e g i e s should be to c o l l e c t data i n key l o c a t i o n s of the flow domain and to i d e n t i f y the s p a t i a l v a r i a t i o n i n h y d r a u l i c c o n d u c t i v i t y . T h i s r e s u l t i s s i m i l a r to that of Smith and Schwartz [19816], The root mean square e r r o r i s a b e t t e r c r i t e r i o n f o r judging the u n c e r t a i n t y of d i s c h a r g e p r e d i c t i o n s than the c o e f f i c i e n t of v a r i a t i o n . The c o e f f i c i e n t of v a r i a t i o n i s a good c r i t e r i o n of accuracy only when the standard d e v i a t i o n gets smaller as the b i a s gets s m a l l e r . For the estimates of pond d i s c h a r g e , the b i a s can be very low with a high standard d e v i a t i o n . Or c o n v e r s e l y , the b i a s can be high with a low standard d e v i a t i o n . For the re f e r e n c e cases c o n s i d e r e d , an i n i t i a l borehole l o c a t e d below the pond tended to be the p r e f e r r e d l o c a t i o n . Sampling i n the region below the pond was p r e f e r r e d more o f t e n than sampling i n other r e g i o n s . For two boreholes, the r e s u l t s suggest that l o c a t i n g both boreholes under the pond i s a b e t t e r s t r a t e g y f o r the SUMMARY / 151 cases t e s t e d than l o c a t i n g one below the pond and one elsewhere. When m u l t i p l e boreholes are d i s t r i b u t e d evenly through the flow domain, i n c r e a s i n g the number of boreholes does not n e c e s s a r i l y l e a d to more c e r t a i n p r e d i c t i o n s of pond d i s c h a r g e or to lower va l u e s of the o b j e c t i v e f u n c t i o n . Whether c o l l e c t i n g 30 samples from 10 boreholes i s p r e f e r a b l e to c o l l e c t i n g 30 samples from 6 boreholes i n s t r a t i f i e d media depends upon the value of the cost c o e f f i c i e n t f o r the l o s s f u n c t i o n and the s p a t i a l v a r i a t i o n i n h y d r a u l i c c o n d u c t i v i t y . In s t r a t i f i e d media, although measurements reduce the u n c e r t a i n t y i n h y d r a u l i c c o n d u c t i v i t y over a l a r g e h o r i z o n t a l d i s t a n c e , c o n s i d e r a b l e u n c e r t a i n t y i n disc h a r g e p r e d i c t i o n s can s t i l l e x i s t even with a r e l a t i v e l y l a r g e number of measurements. When the h y d r a u l i c c o n d u c t i v i t y measurements are spaced evenly through the flow domain, important shallow l a y e r s may be missed that have a l a r g e i n f l u e n c e on pond d i s c h a r g e . RECOMMENDATIONS 1. A l a r g e number of re f e r e n c e cases would be r e q u i r e d to ob t a i n p r o b a b i l i t y d i s t r i b u t i o n s on the p r e d i c t i o n u n c e r t a i n t y f o r a sampling s t r a t e g y and to determine an SUMMARY / 152 optimal sampling s t r a t e g y . T h i s i s because there i s s i g n i f i c a n t v a r i a t i o n i n p r e d i c t i o n u n c e r t a i n t y between re f e r e n c e cases. 2. Even a l a r g e number of measurements of h y d r a u l i c c o n d u c t i v i t y may l e a d to poor p r e d i c t i o n s of d i s c h a r g e . An a l t e r n a t i v e approach to a d e t a i l e d sampling program i s to conduct a pump t e s t . Future s t u d i e s c o u l d compare the c o s t s and p r e d i c t i v e a b i l i t y of these two approaches. 3. For f u t u r e s t u d i e s of t h i s type, the importance of a s t o c h a s t i c f r e e s u r f a c e model should be c o n s i d e r e d . If the water t a b l e i s allowed to f l u c t u a t e in response to v a r i a t i o n s in h y d r a u l i c c o n d u c t i v i t y , p r e d i c t i o n s of d i s c h a r g e might be d i f f e r e n t . These r e s u l t s demonstrate the i n f l u e n c e of s p a t i a l v a r i a t i o n i n h y d r a u l i c c o n d u c t i v i t y on p r e d i c t i o n s of groundwater d i s c h a r g e . A sampling s t r a t e g y may y i e l d a low value of the o b j e c t i v e f u n c t i o n f o r one r e f e r e n c e case but a high v a l u e f o r another. A r e f e r e n c e case can be thought of as analogous to a f i e l d s i t e . Although a sampling s t r a t e g y may have a h i g h p r o b a b i l i t y of l e a d i n g to a good p r e d i c t i o n of d i s c h a r g e , i t can l e a d to poor p r e d i c t i o n s of d i s c h a r g e f o r a p a r t i c u l a r case. BIBLIOGRAPHY Andersson, J . , A. M. Shapiro, and J . Bear, A s t o c h a s t i c model of a f r a c t u r e d rock c o n d i t i o n e d by measured i n f o r m a t i o n , Water Resour. Res., 20(1), 79-88, 1984. Bakr, A. A., L. W. Gelhar, A. L. Gutjahr, and J . R. MacMillan, S t o c h a s t i c a n a l y s i s of s p a t i a l v a r i a b i l i t y i n subsurface flows, 1, Comparison of one- and three-di m e n s i o n a l flows, Water Resour. Res. , 1 4 ( 2 ) , 263-271, 1978. Benjamin, J . R., and C. A. C o r n e l l , P r o b a b i l i t y , St at i s l i c s , and Decision for C i v i l Engineers, McGraw H i l l , N. Y., 1970. Bogardi, I., A. Bardossy, and L. Duckstein, M u l t i c r i t e r i o n network design using g e o s t a t i s t i c s , Water Resour. Res., 2 1 ( 2 ) , 199-208, 1985. C a r r e r a , J . , E. Usunoff, F. Szid a r o v s z k y , A method f o r optimal o b s e r v a t i o n network design f o r groundwater management, /. of Hydrology, 73, 147-163, 1984. C l i f t o n , P. M., and S. P. Neuman, E f f e c t s of k r i g i n g and in v e r s e modeling on c o n d i t i o n a l s i m u l a t i o n of the Avra v a l l e y a q u i f e r i n southern A r i z o n a , Water Resour. Res., 18(4), 1215-1234, 1982. G., S t o c h a s t i c modeling of groundwater flow by u n c o n d i t i o n a l and c o n d i t i o n a l p r o b a b i l i t i e s , 1, C o n d i t i o n a l s i m u l a t i o n and the d i r e c t problem, Water Resour. Res., 18(4), 813-833, 1982. G., S t o c h a s t i c modeling of groundwater flow . by u n c o n d i t i o n a l and c o n d i t i o n a l p r o b a b i l i t i e s : The in v e r s e problem, Water Resour. Res., 27(1), 65-72, 1 985. Davis, J . C , S t a t i s t i c s and Data Analysis in Geology, John Wiley & Sons, N. Y., 1973. Delhomme, J . P., S p a t i a l v a r i a b i l i t y and u n c e r t a i n t y i n groundwater flow parameters: A g e o s t a t i s t i c a l approach, Water Resour. Res., 15(2), 269-280, 1979. Delhomme, J . P., K r i g i n g i n the hydro s c i e n c e s , i n Flow Through Porous Media: Recent Developments, e d i t e d by G. F. Pinder, pp. 99-113, CML p u b l i c a t i o n s , U.K., 1 983. Dagan, Dagan, 153 / 154 D e t t i n g e r , M. D., and J . L. Wilson, F i r s t order a n a l y s i s of u n c e r t a i n t y i n numerical models of groundwater flow, 1, Mathematical development, Water Resour. Res., 77(1), 149-161, 1981. E l - K a d i , A. I., and W. B r u t s a e r t , A p p l i c a b i l i t y of e f f e c t i v e parameters f o r unsteady flow i n nonuniform a q u i f e r s , Water Resour. Res., 21(2), 183-198, 1985. Freeze, R. A., A s t o c h a s t i c - c o n c e p t u a l a n a l y s i s of one-dimensional groundwater flow i n nonuniform homogeneous media, Water Resour. Res., 11(5), 725-741, 1975. Freeze, R. A., A s t o c h a s t i c - c o n c e p t u a l a n a l y s i s of r a i n f a l l - r u n o f f processes on a h i l l s l o p e , Water Resour. Res., 16(2), 391-408, 1980. Freeze, R. A. and J . A. Cherry, Groundwater, P r e n t i c e - H a l l , N. J . , 1979. Freund, J . E., Mathematical Statistics, P r e n t i c e - H a l l , N. J . , 1971. F r i n d , E. 0., and G. B. Matanga, The dual f o r m u l a t i o n of flow f o r contaminant t r a n s p o r t modeling, 1, Review of theory and accuracy a s p e c t s , Water Resour. Res., 21(2), 159-169, 1985. F r i n d , E. 0., G. B. Matanga, and J . A. Cherry, The dual f o r m u l a t i o n of flow f o r contaminant t r a n s p o r t modeling, 2, The Borden a q u i f e r , Water Resour. Res., 21(2), 170-182, 1985. Gelhar, L. W., S t o c h a s t i c subsurface hydrology from theory to a p p l i c a t i o n s , Water Resour. Res., 22(9), 135S-145S, 1986. Hughes, J . P. and D. P. Lettenmaier, Data requirements f o r k r i g i n g : E s t i m a t i o n and network d e s i g n , Water Resour. Res., 17(6), 1641-1650, 1981. J o u r n e l , A. G., and C. J . H u i j b r e g t s , Mining Geostatistics, Academic P r e s s , N. Y., 1978. K i t a n i d i s , P. K., and E. G. Vomvoris, A g e o s t a t i s t i c a l approach to the in v e r s e problem i n groundwater modeling (steady s t a t e ) and one-dimensional s i m u l a t i o n s , Water Resour. Res., 19(3), 677-690, / 155 1983. Loague, K. M., An assessment of r a i n f a l l - r u n o f f modeling methodology, Ph.D. t h e s i s , Dep. of G e o l . S c i . , Univ. of B r i t . Columbia, Vancouver, Canada, 1986. Maddock, T., I I I , Management model as a t o o l f o r studying the worth of data, Water Resour. Res., 9(2), 270-280, 1973. Mantoglou, A. and J . L. Wilson, S i m u l a t i o n of random f i e l d s with the t u r n i n g bands method, Ralph M. Parsons Laboratory, Hydrology and Water Resources Systems, Dept. of Civil Eng., Report No. 264 R81-16,, Mass. I n s t i t u t e of Technology, Cambridge, Mass., 1981. Mantoglou, A. and J . L. Wilson, The t u r n i n g bands method f o r s i m u l a t i o n of random f i e l d s using l i n e g e n e r a t i o n by a s p e c t r a l method, Water Resour. Res., 75(5), 1379-1394, 1982. Massmann, J . and R. A. Freeze, Groundwater contamination from waste management s i t e s : The i n t e r a c t i o n between r i s k - b a s e d e n g i n e e r i n g design and r e g u l a t o r y p o l i c y , 1, Methodology, Water Resour. Res., (in press). Massmann, J . and R. A. Freeze, Groundwater contamination from waste management s i t e s : The i n t e r a c t i o n between r i s k - b a s e d e n g i n e e r i n g design and r e g u l a t o r y p o l i c y , 2, R e s u l t s , Water Resour. Res., (in press). Matheron, G., The i n t r i n s i c random f u n c t i o n s and t h e i r a p p l i c a t i o n s , Advan. Appl. Prob., 5, 439-468, 1973. M e j i a , J . M. and I. Rodriguez-Iturbe, On the s y n t h e s i s of random f i e l d sampling from the spectrum: An a p p l i c a t i o n to the g e n e r a t i o n of h y d r o l o g i c s p a t i a l processes, Water Resour. Res., 70(4), 705-711, 1974. Nakamoto, D. B., F. B. McLaren, and P. J . P h i l l i p s , M u l t i p l e completion monitor w e l l s , Groundwater Monitoring Review, ( S p r i n g ) , p. 50-55, 1986. Neuman, S. P., S t a t i s t i c a l c h a r a c t e r i z a t i o n of a q u i f e r h e t e r o g e n e i t i e s : An overview, G. S. A. Special Paper, 189, 81-102, 1982. Ostendorf, D. W., Modeling contamination of shallow unconfined a q u i f e r s through i n f i l t r a t i o n beds, Water Resour. Res., 22(3), 375-382, 1986. / 156 Rodriguez-Iturbe, I . and J . M. M e j i a , The design of r a i n f a l l networks in time and space, Water Resour. Res., 70(4), 713-728, 1974. Rouhani S., Variance r e d u c t i o n a n a l y s i s , Water Resour. Res., 27(6), 837-846, 1985. Sagar, B., G a l e r k i n f i n i t e element procedure f o r a n a l y s i n g flow through random media, V/ater Resour. Res., J4(6), 1035-1044, 1978. Scheuer, E. M. and D. S. S t o l l e r , On the gene r a t i o n of normal random v e c t o r s , Technometri cs, 4, 278-281, 1 962. Smith, L., S p a t i a l v a r i a b i l i t y of flow parameters i n a s t r a t i f i e d sand, Mathematical Geology, 73(1), 1-21, 1981 . Smith, L., S t o c h a s t i c models of f l u i d flow i n heterogeneous media, ISSS Symposium on Soil Spatial Variability, Las Vegas, Nevada, Nov. 30- Dec. 1, 1984. Smith, L. and R. A. Freeze, S t o c h a s t i c a n a l y s i s of steady s t a t e groundwater flow i n a bounded domain, 1, One-dimensional s i m u l a t i o n s , Water Resour. Res., 75(3), 521-528, 1979a. Smith, L. and R. A. Freeze, S t o c h a s t i c a n a l y s i s of steady . s t a t e groundwater flow i n a bounded domain, 2, Two-dimensional s i m u l a t i o n s , Water Resour. Res., 75(6), 1543-1559, 1979*. Smith, L. and F. W. Schwartz, Mass t r a n s p o r t , 2, A n a l y s i s of u n c e r t a i n t y i n p r e d i c t i o n , Water Resour. Res., 77(2), 351-369, 1981a. Smith, L. and F. W. Schwartz, Mass t r a n s p o r t , 3, Role of h y d r a u l i c c o n d u c t i v i t y data i n p r e d i c t i o n , V/ater Resour. Res., 77(5), 1463-1479, 19816. Swaney, J . M. and E. H. Jackson, Rapid i n f i l t r a t i o n e f f l u e n t c o l o u r removal system at Crestbrook's Skookumchuck p u l p m i l l , B. C. Professional Engineer, (May), p. 21, 1 983. Townley, L. R., Second order e f f e c t s of u n c e r t a i n t r a n s m i s s i v i t i e s on p r e d i c t i o n s of p i e z o m e t r i c heads, Proceedings of the 5th I n t e r n a t i o n a l / 157 Conference on F i n i t e Elements in Water Resources, B u r l i n g t o n , Vermont, p. 251-264, June 1984. Warren, J . E. and H. S. P r i c e , Flow i n heterogeneous porous media, Soc. of Petrol. Eng. J., 7(3), 153-169, 1961. Witherspoon, P. A. and T. N. Narasimhan, G e o t e c h n i c a l e n g i n e e r i n g : I n v e s t i g a t i o n s of a r t i f i c i a l recharge and p o l l u t i o n from l a n d f i l l o p e r a t i o n s , Hydrologic Engineering Center Contract No. DACWO5-71-C-0114, Dept. of C i v i l Eng., Univ. of C a l i f . , Berkeley, 1 972. APPENDIX : TABLE OF NOTATION C / l o s s c o e f f i c i e n t [$T/L2] l o s s f u n c t i o n [$] C o cost of sampling [$] C{ cost of a l t e r n a t i v e treatment [ $ / L 3 ] K h y d r a u l i c c o n d u c t i v i t y , [ L ] / [ T ] Lw width of the pond [L] M decomposed c o v a r i a n c e matrix AC number of Monte C a r l o runs p number of s t o c h a s t i c blocks Q estimate of flow through pond f o r a r e a l i z a t ion Q average estimate of flow through pond f o r a s i m u l a t i o n Qt a c t u a l flow through the pond f o r r e f e r e n c e case r c o r r e l a t i o n c o e f f i c i e n t RCN r e f e r e n c e case number S standard d e v i a t i o n of estimates of flow through the pond from a s i m u l a t i o n Sc standard d e v i a t i o n of l o g K from the set of r e a l i z a t i o n s forming a c o n d i t i o n a l s i m u l a t i o n S™ estimate of standard d e v i a t i o n i n l o g K from a set of borehole measurements 158 / 159 standard d e v i a t i o n of l o g K i n a l l b l o c k s forming the r e f e r e n c e case standard d e v i a t i o n of l o g K from the set of r e a l i z a t i o n s forming an u n c o n d i t i o n a l simulat ion time h o r i z o n [T] c o v a r i a n c e matrix log h y d r a u l i c c o n d u c t i v i t y average value of l o g K from the set of r e a l i z a t i o n s forming a c o n d i t i o n a l s i m u l a t i o n estimates of mean l o g K from set of borehole measurements average value of l o g K i n a l l blocks forming the r e f e r e n c e case average value of l o g K from the set of r e a l i z a t i o n s forming an u n c o n d i t i o n a l s i m u l a t i o n o b j e c t i v e f u n c t i o n [$] i n t e g r a l s c a l e i n d i r e c t i o n x i n t e g r a l s c a l e i n d i r e c t i o n z ensemble mean l o g K a u t o c o r r e l a t i o n ensemble standard d e v i a t i o n l o g K / 160 h y d r a u l i c head, [L]/[T] ii stream f u n c t i o n , [L2]/[T]