S I M U L T A N E O U S I N V E R S I O N O F T H E R M A L A N D H Y D R O G E O L O G I C D A T A by A L L A N D A V I D W O O D B U R Y M . S c , The University of Bri t ish Columbia , 1983 A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S The Department of Geological Sciences We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A A p r i l 1987 © A l l a n David Woodbury, 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 The University of British Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 DE-6(3/8-n To Gunny A B S T R A C T The question that is addressed in this thesis is: can a simultaneous inverse scheme involving thermal and hydrologic data resolve hydrologic model parameters to a better degree than hydrologic data alone? The first chapter sets the framework for this question by first reviewing linear and non-linear inverse problems and then illustrating the advantages of a simultaneous inverse of two different data sets through the use of a simple example. It is the goal of Chapter 2 to examine current methodologies for stating and solving the inverse problem. A review of the maximum likelihood approach is presented, and a construction formalism is adopted by introducing a series of objective functionals (norms) which are minimized to yield a variety of possible models. The inverse is carried out using a modification of a constrained simplex procedure. The algorithm requires no derivative computations and can be used to minimize an arbitrarily complicated non-linear func-t ional , subject to non-linear inequality constraints. The algorithm produces a wide variety of acceptable models clustered about a global minimum, each of which generates data that match observed values. The inverse technique is demonstrated on a series of one and two-dimensional synthetic data sets, and on a hydraulic head data set from Downie Slide, B r i t i s h Columbia , Canada. At this site, four parameters are determined; the free-surface posit ion of the water table and three boundary conditions for the domain. Further s im-ulations using a theoretical data set with assumed properties similar to that of Downie Slide show that with noise free data, and an adequate spacing between points it is possible to interpolate an unbiased estimate of hydraulic head data at all nodes in the equivalent discretized domain. When the inverse technique is applied, the domain's conductivity structure is correctly identified when enough prior log-conductivity information is avail-able. The implications for Downie Slide are that in order to construct anything but a simple hydrogeologic model, accurate field measurements of hydraulic head are required, as well as well-defined estimates of hydraulic conductivity, a better spacing between mea-surements, and adequate knowledge of the boundary conditions. i Chapter 3 is devoted to developing the idea of a joint inversion scheme involving both thermal and hydrologic data. One way of overcoming data l imitations (sparse hydraulic head or few hydraulic conductivity estimates) in an inverse problem is to introduce an inde-pendently collected data set and apply simultaneous or joint inversion. The joint inversion method uses data from a number of different measurements to improve the resolution of parameters which are in common to one or more functional relationships. One such data set is subsurface temperature, which is sensitive to variations in hydraulic conductivity. In Chapter 3, the basic concepts of heat and fluid transfer in porous media wi th emphasis on forced convective effects are reviewed, followed by inversion of theoretical data and a re-investigation of the hydrogeology of Downie Slide, augmented wi th thermal data and a simultaneous inverse. Addi t iona l runs on a heterogenous medium presented in Chapter 2 are carried out. W i t h a good temperature data base, thermal properties can be prop-erly resolved. However, in this stochastic problem the addition of thermal data did not condition .the inverse to a greater degree than accomplished wi th the addition of prior information on log-conductivity. The benefits of including thermal data and applying a joint inversion can be substantial when considering the more realistic problem of uncertain boundary conditions. The simultaneous inverse is also applied to the Downie Slide data set examined in Chapter 2. Unfortunately, wi th a homogeneous hydraulic conductivity, all that can be determined from a hydraulic head inverse are ratios of flux to hydraulic conductivity. B y including thermal data, the value of hydraulic conductivity can be deter-mined at this site. Some of the model parameters (basal heat flux, thermal conductivity, specified head boundaries) are not resolved well by the joint scheme. None-the-less the constructed models do offer valuable insight into the hydrogeology of the field site. The constructed models persistently show a hydraulic conductivity value of about 1 x 10~ 7 m/sec, which is consistent wi th previous estimates of hydraulic conductivity at the site. A further comparison with the inverse results in Chapter 2 show good agreement between the two inverses for the hydraulic properties. ii T A B L E OF CONTENTS Page Abstrac t i Acknowledgements v i i i C H A P T E R 1 - Overview 1.1 Introduction 1 1.2 Review of Linear and Non-linear Inverse Theory 4 1.3 Linear Inverse Problem (Accurate Data) 6 1.4 Non-linear Inverse 10 1.5 Random Search Methods 13 1.6 Inverse Theory In Groundwater Hydrology 14 1.7 Condi t ional Simulations 16 1.8 Aquifer Inverses 17 1.9 Other Inversion Techniques 21 1.10 Overview of Chapters Two and Three 24 References 27 Figures - 3 0 C H A P T E R 2 - Inversion of Hydraul ic Head Data Summary 32 2.1 Introduction 34 2.2 Literature Review 35 2.3 Functional Relationships and Governing Equations 37 2.4 Construct ion 38 2.5 Example Construction Problem 42 2.6 Opt imiza t ion Technique 44 2.7 Stopping Cr i te r ia 48 iii 2.8 Computat ion of Covariances 50 2.9 Selection of A o p t 51 2.10 Summary of Algor i thm 52 2.11 Inversion of Theoretical Data 54 2.12 Analysis of Field Data 57 2.13 Discussion and Conclusions 69 Appendix 2A 72 Appendix 2B 74 Notat ion 76 References 79 Tables 83 Figures 93 C H A P T E R 3 - Hydraulic and Thermal Inversion Summary I l l 3.1 Introduction 113 3.2 Literature Review 116 3.3 Functional Relationships and Governing Equations 117 3.4 Construction 121 3.5 Inversion Technique 122 3.6 Summary of Algor i thm 123 3.7 Inversion of Theoretical Data 124 3.8 Analysis of Field Data 136 3.9 Discussion and Conculsions 146 Notat ion 149 References 152 Tables 154 Figures 160 C H A P T E R 4 - Conclusions Summary 177 iv LIST OF TABLES Table Page 2.1. Example Construction Problem 1 83 2.2. Example Construction Problem 2 84 2.3. Example Construction Problem 3 85 2.4. One-dimensional, steady state inverse 86 2.5. Comparison between Carrera and Neuman 87 2.6. Downie Slide hydraulic head data set 88 2.7. Comparison of kriging results 89 2.8. Input parameter ranges for inverse 90 2.9. Comparison of results for Downie Slide inverse 91 2.10. Synthetic hydraulic head data set 92 3.1. Synthetic hydraulic head and temperature data set 154 3.2. Sparse synthetic hydraulic head and temperature data set 155 3.3. Synthetic hydraulic head and temperature data set 156 3.4. Comparison of actual and inverse results 157 3.5. Input parameters and ranges 158 3.6. Downie Slide inverse results 159 v LIST OF ILLUSTRATIONS Figure Page 1.1. Graph of a straight line 30 1.2. Graph of a curve 31 2.1. Graph of a curve and four data points 93 2.2. One dimensional groundwater flow system 94 2.3. Objective function J as a function of model parameters 95 2.4. Log-conductivity realizations 96 2.5. Example two-dimensional inverse problem 97 2.6. Regional groundwater flow system 98 2.7. Site map of Downie Slide 99 2.8. Fini te Element mesh 100 2.9. Variogram for vertically oriented points 101 2.10. Var iogram for scattered data 102 2.11. Kr iged theoretical hydraulic head data 103 2.12. Downie Slide inverse problem 104 2.13. Es t imat ion errors in hydraulic head data 105 2.14. Li inverse hydraulic heads 106 2.15. L 2 inverse standard deviations 107 . 2.16. Log-conductivity field L2 inverse, prior information 108 2.17. Ac tua l log-conductivity field 109 2.18. Smallest deviatoric log-conductivity field 110 3.1. Thermal effects of groundwater flow 160 3.2. Fini te Element mesh 161 3.3. Objective function J as a function of model parameters 162 3.4. Four temperature fields 163 3.5. Ac tua l model for snythetic example 164 3.6. Hydraul ic head field 165 vi 3.7. Temperature field for shallow aquifer 166 3.8. Locat ion predicted by inverse (heads) 167 3.9. Locat ion predicted by inverse (temps) 168 3.10. Groundwater flow domain 169 3.11. Locat ion of aquifer predicted by inverse 170 3.12. Locat ion of aquifer predicted by inverse 171 3.13. Log-conductivity structure 172 3.14. Downie Slide site map 173 3.15. Thermal profiles 174 3.16. Downie Slide inverse 175 3.17. Temperature field pedicted from inverse 176 vii A C K N O W L E D G E M E N T S Words alone cannot describe my appreciation to the individuals and organizations who contributed to this work. I would like to thank my thesis supervisor, Leslie Smith, for his tremendous support and reviews of my work. M y other advisors were Scott Dunbar, A l Freeze, Mike Novak, and Doug Oldenburg. I would like to thank all these people, and in particular Scott Dunbar, for the many interactions that we have had. Randy Bourne at B . C . Hydro is also deeply thanked for his assistance throughout the project. Tom Nicol and J i m Glosl i at the U . B . C . computing centre helped considerably in finding bugs in my computer codes. Bo Chandra unselfishly gave time to help construct a workable thermal probe. Various pieces of equipment were supplied by Westbay Instruments L t d . and Nevin Sadler Brown Goodbrand Inc. Addi t iona l field assistance was provided by Gregg Dallas, Rene Therien, and McElhanney Engineering. This research was funded by grants from the Natural Sciences and Engineering Research Counci l of Canada ( N S E R C ) . Computa-tions were carried out on a F P S 1 6 4 / M A X array processor supported by an N S E R C major installation grant to the University of Br i t i sh Columbia . I would finally like to express my appreciation to the management of Br i t i sh Columbia Hydro and Power Authority for permission to publish the material on Downie Slide. viii C H A P T E R 1 Overview 1.1 Introduction The difficulties associated wi th making direct measurements of all the hydrologic pa-rameters associated wi th physically based mathematical models are well known. Equal ly well known are the difficulties in trying to adjust model parameters unti l hydraulic head output at selected points match observed values. Cal ibrat ion can be a time consuming, expensive, and frustrating experience. Quite frequently, questions are raised as to whether or not the final model parameters represent opt imum values, and to what degree the model is unique. A classic example of the 'inverse problem' or ' inversion' is the determination of the transmissivity (T) and storativity (5) of an aquifer from hydraulic head measurements taken in observation wells adjacent to a pumping well . In this case, hydraulic heads constitute a set of observed data, and the idealized case of an infinite, horizontal aquifer (depending on the boundary value problem) with constant T and S values constitutes the 'model ' . In this thesis, the term 'model ' is defined in a different sense than what the reader is probably accustomed to, and is defined as a configuration of material properties. In some cases, boundary conditions are also treated as components of the model. A n estimate of the model produces a set of calculated hydraulic heads at the observation wells. The inevitable differences between the observed and calculated hydraulic heads at these points can be attr ibuted to both random and systematic measurement errors, and blunders. The sum total of all these errors is usually referred to as 'noise'. A n oversimplification of the physics incorporated in the boundary value problem which forms a basic component of 1 1.1 Introduction 2 the inverse can also lead to another component of noise. Usually, estimates of varying reliabil i ty of T and S in the aquifer are available before a well test. The ranges on these values constitute model constraints. In the usual approach, T and S are found by a tr ial and error curve matching procedure in which drawdown curves are adjusted to match a set of type curves; the closeness of the fit determines the 'resolution' of the answer. Unfortunately, the calculated values of T and S are sensitive to the exact form of the type curve. Different type curves result from different boundary value problems, for example; those that incorporate leakage from confining layers. In other words, a slightly different type curve wi l l yield values of T and S that may vary by orders of magnitude. This non-uniqueness illustrates a major problem wi th inverse theory, namely that errors in physical assumptions can lead to vast differences in models which are constructed. Al though this overview is not concerned wi th the inversion of well test data, the concepts introduced here serve as an analogy to the generalized inverse theory that wi l l be presented in the next sections. This chapter sets the framework for simultaneous inversion of hydraulic head and temperature data by reviewing the linear and non-linear inverse problem and illustrat-ing the advantages of a simultaneous inverse of two different data sets through the use of a simple example. Al though inversion theory has been the focus of many papers in the geophysical and hydrogeologic literature (for example, gravity, seismic, resistivity, magnetotelluric, magnetic, and conductive heat flow), little work has been reported on simultaneous inversion. Most of the research in this area has been limited to combined resistivity/electromagnetic problems (Petrick et. a l , 1977; Ward et. al , 1976; Glenn and Ward , 1976). In groundwater hydrology, inversion has been mainly applied to finding best-fitting hydraulic conductivity distributions in aquifers given hydraulic head data. A recent review of this inverse problem is given by Yeh (1986). It has been shown (for example, Smi th and Chapman, 1983) that groundwater can play an important role in disturbing the near-surface thermal regime. These authors have 1.1 Introduction 3 demonstrated the sensitivity of the thermal field to: • Variations in the magnitude of permeability and anisotropic ratios. • Length /depth ratios of the flow system. • Thermal conductivity changes. • Variable boundary conditions. • The existence and locations of aquifers. In addit ion, Smi th and Chapman (1983) detailed the sensitivity of an advective threshold. A t low permeability the thermal regime of a basin is purely conductive, because groundwa-ter velocities are too low to effect redistribution of heat. The transition from a conductive to an advectively-disturbed thermal regime is sharp, occurring over one order of magnitude in permeability for the basins they considered. Ult imately, at high fluid velocities, a basin can be advectively dominated wi th a near isothermal temperature field. The governing equation of energy transport in porous media is coupled to the fluid flow equation through the fluid's specific discharge. This physical coupling between the two equations produces a potential benefit. The question that is addressed in this thesis is: w i th in certain ranges of advective disturbance, can a simultaneous inverse scheme involving thermal and hydrologic data resolve hydrologic model parameters to a better degree than hydrologic data alone? If so, the additional data required by the joint inversion scheme would be simple to obtain. Temperature measurements would be taken in piezometers or boreholes along wi th hydraulic head measurements. Temperature measuring equipment is also simple to construct and use in the field. Note that while the number of hydraulic head measurements is l imited in an individual piezometer, temperatures can be profiled continuously along the length of a borehole or piezometer standpipe. This feature greatly increases the amount of temperature data compared to hydraulic heads that can be gath-ered at most sites. The most important l imitat ion remains; that is, l imits in the number of boreholes. It is the goal of this thesis to demonstrate that the joint inverse scheme can lead to substantial improvements in model resolution, particularly in cases where limita-1.2 Review of Linear and Non-linear Inverse Theory 4 tions in hydraulic head data sets, such as incompleteness or inaccuracy, are pronounced. Only steady-state hydraulic head and temperature fields wi l l be considered in this study. Th i s approximation is appropriate where yearly fluctuations in the water table are small compared to the depth of the flow system and where seasonal temperature variations effect only the near-surface temperature field. Th i s chapter reviews the following topics: 1. Basic theory in linear and non-linear inversion. 2. Applicat ions of inverse theory in groundwater hydrology. 3. The work to be carried out in Chapters 2 and 3 in this thesis. 1.2 Review of Linear and Non-linear Inverse Theory The foundations of inverse theory can be found in geophysics, applied mathematics, linear algebra, and electrical engineering in papers and textbooks published in the 1960's and early 1970's. However, the theory developed in geophysics and linear algebra has the most relevancy and similarity to the groundwater problem. As such, inverse theory in this thesis w i l l be developed from a geophysical perspective. Parallel developments in theory as applied to groundwater hydrology are discussed in the next section. A review of the theory of inversion can be found in Wiggins (1972). The primary source of the theoretical development was established by Backus and Gilbert (1967,1968), in their delineation of earth structure from seismic data. In some of the early work in inversion, Inman et al . (1973) dealt wi th the problem of direct interpretation of resistivity curves, and Glenn et al . (1973) presented a discussion of inversion and applied the theory to magnetic dipole data. Oldenburg (1984) provided a good summary and background to the linear problem and Menke (1984) provided a much needed textbook on the subject at the introductory level. Neuman (1973) classified methods of parameter estimation as direct or indirect de-pending on whether they are obtained from a direct 'back calculation' procedure (for 1.2 Review of Linear and Non-linear Inverse Theory 5 example, when all hydraulic heads are known in the flow field), or an indirect calculation resulting from some least-squares optimization scheme. It is the goal of many mathemati-cians working in this field to find a direct solution to the inverse problem. This exercise is useful, part icularly in demonstrating that an inverse is unique when the observed data are perfect, and to aid in conceptual analyses of simple problems. In groundwater hydrology such analytic solutions are not possible for most practical situations, due to the nature of the governing equations and the boundary value problems solved. As such, analytic inverses or direct methods are not discussed in this thesis. However, some numerical -direct methods (Frind and Pinder, 1973; Yeh et al . , 1983) wi l l be discussed in the review of groundwater literature. In order to solve an inverse problem a formalism or philosophy about the problem needs to be developed. For instance, specific goals about what is to be achieved from the inverse need to be stated. These inverse objectives can be generalized into three groups (after Oldenburg, 1984): 1 Construction: Generate a model which reproduces the data. In other words, generate a hydrogeologic structure ( hydraulic conductivities or other parameters) such that the hydraulic head field computed from this structure matches the observations to wi thin a priori standard errors. These standard errors reflect uncertainties in the data set. 2 Appraisal: If the inverse leads to nonunique models, this approach seeks to find unique averages of al l possible models. For instance, this approach would combine several hydraulic conductivity zones for a unique average value. 3 Inference: F i n d absolute upper and lower bounds of the model. When the ability of the data to resolve model parameters is poor, or when the problem is non-unique, this method attempts to find realistic upper and lower limits to the model. It is the primary focus of this thesis to consider only construction methods. In the next few sections the forward problem wi l l be defined, followed by discussions of the linear and non-linear inverse problems. In the following discussions of inverse theory, the assumption 1.3 Linear Inverse Problem (Accurate Data) 6 _ , * • is made that the set of observed data is obtained wi th total accuracy, and that there is a perfect correspondance between physical theory and observation (i.e. no "model" noise). These approximations are made for the sake of clarity, and later in this chapter and in Chapter 2, these restrictions are relaxed in dealing wi th realistic data sets. 1.3 Linear Inverse Problem (Accurate Data) The forward problem is solved by constructing a functional relationship to predict physical data for a given set of input parameters to a physical model. That is, construct a functional relationship: d = *{x,p) (1.1) where d is an M dimensional vector of computed values of some physical quantity, and 5 is a function of the model parameters contained in the vector p. Spatial coordinates are contained in vector x. The inverse problem can be stated as: given an incomplete and inaccurate data set, extract information about the model. ' Information' in this sense refers to either determin-ing any model which reproduces the data, a unique average of the model, or upper and lower bounds to the model. In many applications of inverse theory in geophysics, equation ( l . l ) takes the form of a Fredholm integral equation of the first k ind and the unknown model p is a continuous function of the spatial coordinates. Oldenburg (1984) summa-rized inverse methods and algorithms suited to this type of problem. However, in many instances the continuous model must be replaced by a finite set of parameters. The model p is parti t ioned into segments (i.e. finite elements or some sort of layering ) and more data are collected than unknown parameters sought. This 'parametric ' approach leads to an overdetermined system of equations and a unique solution can be found by least-squares. Another approach must be taken when the number of unknown parameters exceeds the number of observation points. This situation is known as an over-parameterized problem and leads to an underdetermined system of equations. Undertermined systems lead to 1.3 Linear Inverse Problem (Accurate Data) 7 non-unique solutions for model parameters. Groundwater flow domains which are mod-eled in cross section could lead to this situation because of the large spatial variability in hydraulic conductivity that could be possible and that usually few hydraulic head data are available. In contrast, when identifying hydraulic conductivity values for an areal two-dimensional aquifer, an overdetermined system of equations would usually result. Equat ion (1.1) may be linear; for instance, the gravity effect over a sphere: d = Ap (1.2) where A is an M x N matrix of coefficients that represents the physics of the problem and the location of the observation points, M is the number of data points, and p is the model (in this case-unknown density of the sphere), of dimension N, where N is number of model parameters. It is useful at this stage to define the singular value decomposition (SVD) of A. The S V D provides a computationally efficient (and stable) alternative to matrix inversion. A can be writ ten as: A = USVT (1.3) where U is an M x M orthogonal matr ix, V is a N x N orthogonal matrix and S is an M x N diagonal matr ix of real numbers, known as the singular values of A. These values are usually arranged in descending order (Forsythe, et. al , 1977). The reader may recall that the columns of an orthogonal matr ix have unit length and are orthogonal to one another (orthonormal). If the columns of a matrix are orthogonal, then the dot product of any two columns is equal to zero. The result of these features leads to useful property of a matr ix as a whole; namely, the transpose of a matrix is equal to its inverse. For example, uT = u-\ A solution to equation (1.2) is found by minimizing the difference between d and Ap. This difference is defined by means of various norms, which are scalar measures of length. The most common measure of the length of a vector is its Euclidean or Li length defined 1.3 Linear Inverse Problem (Accurate Data) 8 as: M N l l L 3 = ( £ r f ? ) 1 / 2 (1-4) 1=1 Other measures of length are possible, for instance, the L\ norm: M 1=1 If the problem is overdetermined a solution is usually sought that minimizes the squared Li difference between the observed and calculated data ( in vector notation ): min J = (d — Ap)T(d - Ap) (1.6) where J is referred to as an objective function. This solution is commonly called the least-squares solution to (1.2) and can be wri t ten as: pest= (ATA)~lATd =VS~lUTd (1.7) where pest is the estimated model from the inverse. In the case where there is less data than parameters (i.e. M < N), the problem is underdetermined and a unique solution of the least squares problem (1.7) is not possible. A minimum length solution can be obtained by solving the constrained problem: min J = p p subject to: d = Ap (1.8) for which the solution is: p0 = AT[AAT}-1d = VqS~lUjd (1.9) and the q subscript indicates that some of the singular values of the matrix A have been discarded. The problem is then said to have been 'reparameterized'. Note that the rank of A in the underdetermined case can be at most equal to M. Under this scheme, high 1.3 Linear Inverse Problem (Accurate Data) 9 absolute values of model parameters are penalized, such that minimizat ion of the norm finds the lowest absolute model values (within the constraints) that satisfy the data. The general solution of an underdetermined or rank deficient problem is non-unique, and can be represented by: P = Po + P± (1.10) where p is the general solution to the inverse problem, p0 represents the min imum length solution, and p1- is an arbitrary vector of null-space components (p1- = V0a), where V0 is an (N x M — N) matrix of orthonormal vectors associated wi th the zero singular values of A , and a is a ( M — N x 1) vector of arbitrary constants. The inner-product of p 0 and p1- satisfies the relation: ( p o , p x ) = 0 (1.11) The S V D provides a useful general formalism for the solution of linear problems. In addition to the solution to (1.2), a resolution and information density matr ix can be defined. A resolution matr ix can be defined as: R = VqVqT (1.12) This matr ix should be an identity matrix, if the resolution of the parameters is perfect. Resolution can be viewed as how well an inverse can identify or separate physically adjacent model parameters. The closeness of the rows in R to exhibit unity on the diagonals and zeros off the diagonals is taken as a measure of this resolution. As might be expected, there is a classic trade-off between resolution and inverse-model variances. As the resolution of the model parameters increases, then so too does the uncertainty of the estimate (Menke, 1984). A n information density matrix can be defined as: S = UqU? (1.13) This matrix gives an indication of how well the observed data reproduces the actual data values. If the observed data and the inverse reproduce the actual data well , then the 1.4 Non-linear Inverse diagonals of this matrix wi l l be close to a value of one. Departures from diagonal dominance wi l l occur when there is a strong interdependence of information among the data or where there are a number of small derivatives for any parameter at certain data points. These two matrices can be used to provide criteria for the design of data collection programs or in the evaluation of the information from an inverse. 1.4 Non-linear Inverse Accurate Data The governing equations in groundwater hydrology generally yield non-linear func-t ional transformations that relate hydraulic head at a point with an integral expression involving hydraulic conductivity (for instance, a functional expression from a finite element equation). If the governing equations can be made linear, then techniques outlined in O l d -enburg (1984) for a linear inverse can be used. In particular, if the governing equations can be decomposed and writ ten in terms of Green's functions, then the inverse problem becomes relatively straight forward. For a non-linear problem there is no general theory that allows for a complete solution. Usually, a parametric approach is followed, in which a quasi-linear system is obtained by expanding (1.1) about some point p° in parameter space and neglecting higher order terms (see Oldenburg, 1979). For instance, d(p,x): = d(p°,x),- + f > - P j ) ^ ^ i ^ ° + (1.14) 3 = 1 ° P j and where py is the j th unknown parameter, d{p,x)l is an observed value of d t , i = 1,2,.., M is the number of data points. If there are N unknown parameters and if (1.14) is rearranged and writ ten for M data points then a matrix equation results: Ad = AAp ' (1.15) where A is a M x N Jacobian or 'sensitivity' matrix, Ad{ = d(p,x)* — d ( p ° , x ) , , and Apj = (pj - py). Once A p is solved for, it can be added to the starting guess, p ° , to arrive 1.4 Non-linear Inverse at an updated value of the model, p 1 . Al though a linear system of equations exists because higher order terms are eliminated, several iterations of this procedure may be needed to converge to obtain a solution. In order to solve for A p , the difference between Ad and AAp must be minimized in some sense. Wiggins (1972) casts the problem in terms of a generalized inverse that is val id if the matr ix system (115) is over determined, underdetermined, or rank deficient. A solution to the following is sought: Ap=HAd (1.16) where H = {ATA\-XAT or AT[AAT\-X (1.17) The solution is found by S V D and mult ipl icat ion: A P n + l = Hq(n)Ad (1.18) where n is an iteration level. The new parameter estimate pn+\ which minimizes the least squares criterion is: Pn+l = Pn + Apn+l (1.19) and p° is used as a starting value. Inaccurate Data with Prior Information The non-linear inverse problem can be extended to the case where prior statistical information on the model is available and where the data are noisy. The observed data are now distinguished from their ' true' counterpart, and are denoted as d". If prior statistics are available on d* (say known measurement or interpolation errors) then a prior covariance matr ix of the observed data can be denned as cov d*. Similarly, prior information on the model may also be available in the form of a covariance matrix cov p", and estimates or measurements of the model parameters, p*. These two a priori covariance matrices can be combined into a larger partitioned matr ix , cov z. In this section, a generalized 1.4 Non-linear Inverse inverse under the L2 norm is found that incorporates this prior statistical information. T w o vectors are defined ( adapted from Tarantola and Valette, 1982; Menke, 1984): z* = (d*,p*)T (1.20) and z={d,p)T (1.21) where p* is a prior estimate of the model parameters. Min imiz ing : $ = \ z - z'\T{cov z\-l\z - z*\ (1.22) subject to exact theory (i.e. no model noise): f(z) = d - S ( z , p ) = 0 (1.23) leads to the equation (Menke, 1984): Pen%=P'+Gn9(\d'-^(x,penst)} + Gn[Penst-p'}) (1-24) Here n is an iteration index, and G~' = ({cov p * ] - 1 + GTn[cov d*\-lGZ)-lGl[c<» d*]"1 (1-25) where is evaluated at p ^ s t , assuming that the measured hydraulic head and model parameters are uncorrelated. This assumption is reasonable since hydraulic conductivity and hydraulic head are measured in fundamentally different ways. G measures the sensitivity of values at data points to changes in model parameters. Individual terms in this matrix can be calculated by perturbing the model about some point, and approximating any G,-y value by 1.5 R a n d o m Search Methods finite difference. Neuman (1980) adopted an adjoint theory approach developed in other disciplines to compute these sensitivity terms. If one uses (1.24) and if the init ial guess is close enough to the global minimum, the successive approximations wi l l converge. Because the distribution for pest is in general non-Gaussian, the covariance of the model parameters may be difficult to interpret, especially in terms of confidence intervals. This matr ix, (cov p e s t ) , can be approximated as: (Menke, 1984): [cov pest] = G~9[cov d'}G-'T+ [I - Rn}[cov p*][I - Rn]T (1.27) where Rn = G~gGn, I is the identity matr ix, and the last n is used. 1.5 Random Search Methods In the last section, a quasi-linear system of equations was solved iteratively for the nonlinear inverse. Methods have been proposed to solve the inverse in its original, nonlinear form. Wiggins (1969,1972) outlined a Monte-Carlo approach to the inversion of seismic data. This type of procedure is a formalized trial and error approach for finding models that are consistent wi th the observations. The procedure consists of generating random models wi th in constraint limits and testing the calculated values at the observation points. A l l models in which the calculated values are within some multiple of the standard deviation of the observed values are saved. The result is a scatter diagram of the set of passing models. Wiggins (1972) shows that the Monte Car lo approach, wi th an adequate sampling of random models which are solutions to the inverse problem, produces all the quantities needed to satisfy the four general goals of an inverse that he poses: an average solution, the standard error of the parameter estimates, the resolving power in parameter space, and the information distribution among the observations (see section 1.3 for definitions). The most powerful and general techniques to solve the non-linear inverse were devel-oped in the field of chemical engineering. Some of these techniques are search procedures that involve repetitive forward computations. The advantage of these methods is that 1.6 Inverse Theory in Groundwater Hydrology 14 they are the most general for minimizing (or maximizing) an arbitrary objective function; however, the drawback is in computational effort. The inverse problem is posed as a non-linear optimizat ion problem wi th linear or non-linear constraints. Dunbar (1977) used a constrained simplex procedure to solve a non-linear problem in identifing the geometry of a fault, and Silva and Hohmann(1983) used a controlled random search procedure to optimize a non-linear function arising in the interpretation of a magnetic anomaly. These algorithms easily handle constraints, do not require partial derivatives of the functional wi th respect to the parameters, and can be employed to minimize any arbitrarily compli-cated non-linear functional which may contain more than one system of equations. Besides finding the global minimum, these algorithms can be adapted to produce a sampling of the parameter space around the global min imum. This sampling can produce the parameter and estimated data covariances matrices. 1.6 Inverse Theory in Groundwater Hydrology The first application of inverse theory in groundwater hydrology can be found in the work of Nelson (1968), Emsellem and de Mars i ly (1971), and Fr ind and Pinder (1973). As the majority of inversion work has been directed towards aquifer problems, and therefore, most case histories are carried out over areal two-dimensional space. In Emsellem and de Mars i ly ' s paper, the authors used an automated approach in which the standard matrix equation Ah = / in a numerical approximation of the fluid flow equation is rewritten in terms of T being the unknown. Here, A is defined as a global stiffness matr ix , h is a vector of hydraulic heads at nodal points, and / is a vector of imposed boundary conditions. If the values of hydraulic head are known at all nodal points in the discretized flow domain then this formulation is val id . The inverse problem is then linear in T , but potentially non-unique. In order to properly pose the problem mathematically, a boundary condition must be known along a line T cutting every streamline of the flow in aquifer region, or the value of T must be known at one point on each streamline. Either 1.6 Inverse Theory in Groundwater Hydrology of these constraints requires a substantial amount of data. Fr ind and Pinder (1973) used a similar mathematical framework as Emsellem and de Mars i ly except the authors used a finite element solution and pose the problem in a different way. They argued that any pumping well is an example of knowing the value of T along every streamline. Thus, they were able to pose the problem in a correct manner and achieve an inverse. The authors presented a sensitivity example in their paper to show that the computed value of hydraulic head in the aquifer is insensitive to a wide range of T values. However, their solution only involved one unknown hydraulic conductivity parameter, i.e., a homogeneous case. Chang and Yeh (1976) proposed an optimization scheme for the solution of the inverse problem. They viewed parameter identification in terms of an optimization problem which can be solved by quasi-linear programming, and applied the method to a transient problem. The inverse method they proposed identified 84 parameters with relative ease. Their analysis began wi th the transient flow equation: where Q is a source/sink term. The objective is to find the best (optimal) values of system parameters such that a certain criterion is optimized while all other system contraints are satisfied. The constraints in this case are the governing equation and the upper and lower bounds on the parameters. The least squares criterion involves minimizing the objective function: where k is the time level, and h* is the set of observed water levels (assumed noise free), at nodal locations The solution algorithm involves linearizing the original partial differential equation, and uses the finite difference method to solve the final equation, which is wri t ten in terms of an unknown T. Standard quadratic programming is used to obtain a solution. Note that quadratic programming refers to a least squares solution that v-(rvfc) = s— + Q (1.28) (1.29) 1.7 Conditional Simulations involves prior constraints imposed on the answer. The process stops when J falls below a specified value. The only disadvantages noted were the large storage requirement noted for the quadratic programming scheme, and the l imitat ion that only values of T are solved for in the inverse. 1.7 Conditional Simulations In the late 1970's and early 1980's, many hydrologic inverse studies were initiated as a result of Freeze's (1975) and others' work in the application of stochastic techniques in modeling fluid flow in heterogeneous media. In a review paper Neuman (1982) stated that in order to use stochastic approaches a tremendous amount of input data is required. The model is viewed as a random function wi th a mean, variance, and autocovariance. Ear ly stochastic simulations were not conditioned on any measured transmissivity values that may be obtained from an aquifer. Delhomme (1979) developed a stochastic algo-r i t h m that is capable of taking measured transmissivity data into account. It is called a condit ional simulation, and is based on the theory of kriging developed in geostatistics. Gambola t i and Volpi (1979) reviewed kriging in a deterministic sense and show its relation to other interpolative schemes. The advantage of kriging over other interpolation schemes is that estimates of variances in interpolated values can be obtained. In Delhomme's work, however, only transmissivity data were conditioned. The procedure is outlined below: 1 K r i g i n g is used to interpolate log values of transmissivities (as areal or point quantities) z(x) to all locations over the entire aquifer, 2 A realization s(x) having the same statisitics as the kriged transmissivity field is gener-ated wi th a turning bands procedure. This realization may have different values than the original data set at the observation points. 3 Each realization is made consistent wi th the observations. The end result is a condi-tioned realization, Z(x). 1.8 Aquifer Inverses The final stage in the procedure is the computation of hydraulic heads based on each Z[x) realization. The standard deviations in hydraulic head can be computed using a Monte-Carlo technique. Delhomme, on the basis of 50 conditioned runs and 30 log T data points, verified that the kriged map was equal to the mean T map; however, the hydraulic heads did not match the observed heads in some portions of the aquifer. Thus, the available log T data were consistent wi th the values at observation points, but this was not a sufficient cri teria to ensure that a good fit in hydraulic heads would be generated. 1.8 Aquifer Inverses Neuman and Yakowitz (1979) proposed that hydraulic head measurements and trans-missivity measurements should be incorporated in conditional simulations by way of an inverse procedure. In their method, both T and h values are treated as random variables. K r i g i n g is not necessary in the inverse procedure, provided that other estimates of the co-variances of h and T are available. The procedure starts with the standard finite element equation for the steady state problem. A{T)h = Q (1.30) where A is a global stiffness matrix, h is a vector of hydraulic heads at nodal points, and Q is a vector of imposed fluxes. Note that A is a linear function of transmissivity. If h" ,Q~ and T* are measured or interpolated quantities; then, h' = h + n (1.31a) T' = T + v (1.316) Q* = Q + $ ( i . 3 i c ) The error terms n, u and c are assumed to be normally distributed. Substitution of (1.31a) and (1.31c) into the finite element equation results in : A{T)h* = Q* + A{T)n (1.32) 1.8 Aquifer Inverses Here the error in Q* is assumed small in comparison to t j . Rearranging: h~ = {A(T)}-1Q* +r, (1.33) This equation can be viewed as a non-linear regression of h against Q' in terms of model parameters T. The actual estimation of T is accomplished by minimizing the least-squares criteria: J(T) = {h* - h)TV-\h* -h) + \[{T* - TfVf^T* - T)\ (1.34) where and V? are a priori covariance matrices for hydraulic head and transmissivity, h = A~1Q", and T* is a prior estimate of T. The authors presented methods to determine the opt imum value of A, which is a ratio that determines the 'weight' of the transmissivity versus hydraulic head measurements. The system of equations was solved by linearization and a variation of Newton's method. The authors used a numerical approach to compute the sensitivity matr ix (see equation 26). They concluded the paper wi th a theoretical example to demonstrate the use of the theory. This paper was followed by two others in a series. The second paper (Neuman et. a l , 1980), was a case study to illustrate the method proposed in paper 1. Field data (consisting of steady-state hydraulic head and log-conductivity data) from the Cortaro Basin in Ar izona was used in the model study. A t the time, the authors did not have a satisfactory kriging program, so kriging was not used in the study. Instead, the authors estimated nodal values of h" and T~ from interpolated maps. The estimated T"s are shown to be in good agreement with those estimated by a t r ia l and error (calibration) approach. A finite element transient model was then used to successfully reproduce 25 years of pumping in the basin. In the th i rd paper in the series, Neuman (1980) proposed an improved solution tech-nique to the method used in papers 1 and 2. This method is based on an adjoint-state method which allows for the rapid computation of the gradient of the error criterion (see equation 1.34). Neuman also chose to work wi th log T values, which he concluded to be more appropriate. 1.8 Aquifer Inverses Clif ton and Neuman (1982) tested the method in conditional simulations. The authors treated transmissivity as a random variable whose spatial fluctuations are a realization of a stochastic process. Simulations of an aquifer in the A v r a Valley of Ar izona were carried out in three stages. Fi rs t , estimates of head, using an average value of T were computed in an unconditional simulation. This simulation d id not incorporate the spatial relationship of T. The computed hydraulic heads have a relatively large variance. Second, the variance can be reduced by kriging and conditioning the log transmissivity data. The variances in hydraulic head by condit ioning were reduced by a factor of 3.2. T h i r d , the lowest level of uncertainty was found when conditional simulations were tied to the inverse. The result was a reduction in variance of hydraulic head of 46.0. This behaviour points out that the conditioning effect of the inverse is much greater than that based solely on kriging. Further work in this area was carried out by Carrera and Neuman (l986a,b,c). These authors reformulated the inverse in a maximum-likelihood framework for transient simu-lations. In addition, the authors treated uncertainty in boundary conditions in the scheme and applied the technique to a series of synthetic and field examples. In another series of papers Cooley (1977, 1979, 1982, 1983) presented a slightly differ-ent inversion technique to that of Neuman's. Cooley's approach is similar to the param-eterization used in geophysics. Cooley (1977) examined the aquifer inverse as a classical non-linear regression problem wi th the solution of the flow equation forming the regression equation and the hydraulic quantities such as T , S, sources/sinks and boundary conditions as unknown model parameters. In this respect, Cooley's approach was more general than that of Neuman's early work (1977) in that Neuman only considered T as the unknown parameter to be determined. Cooley chose to linearize the governing finite element equa-t ion and use a Gauss-Newton method. The sensitivity matrix (1-26) was developed by a numerical scheme. Cooley (1977) made the following conclusions: 1. Convergence of the scheme is related to the number of unknown parameters. 2. I l l conditioning of the sensitivity matrix results from failure to implici t ly satisfy the 1.8 Aquifer Inverses Cauchy condit ion in the boundary value problem; i.e. failure to have a unique function of T to cross all streamlines. 3. Poor parameter estimates result when more parameters are sought than possible for a given problem. Cooley (1979) applied his algorithm to case studies in two field areas: Truckee Mead-ows, Nevada, and a cross-section in Hu la Bas in , Israel (both steady-state situations). Statist ical techniques were used to estimate the degree of non-linearity in the models, the goodness of fit to the observed field data, and the reliability of predictions made with the models. It was found that the standard errors for parameters in both field sites (after the inversion process) were large. This factor indicates that the solutions may be non-unique, although the model fit was good in the case of Truckee Meadows and fair for the H u l a Bas in . In spite of the indicated non-uniqueness, predicted drawdown distributions for a gravel pit operation in the Truckee Meadows were not statistically different than the observed values. In Cooley's papers in 1977 and 1979, only hydraulic head data was used in the in-verse. However, if prior measurements are available, i.e. T*, etc., then there should be a way of incorporating these measurements and their reliability into the regression model. Cooley (1982) presented such a method, and defined two levels of prior information; prior information of known statistics and prior information of unknown reliability (informed guesses, estimations, etc.). Cooley's treatment was different from that of Neuman (1980) in that prior information from many parameters, not just T , are included; and a different technique of determining the covariance structure of the model was adopted. The ideal case would be where prior information of the first level is available. Cooley concluded that this is not likely to be the case. Methods for obtaining extensive unbiased information for a l l parameters do not exist at present. In many instances only a best guess may exist for any parameter. However, adding any prior information has the effect of stabilizing the inversion process. In a later paper (Cooley, 1983) his inverse approach was applied 1.9 Other Inversion Techniques 21 to several sites, including Truckee Meadows, Nevada. Results for this area showed that the inclusion of prior information did not change the parameter estimates from that of ordinary regression ( Cooley, 1979), however, analysis showed that more reliability could be placed in the estimates. 1.9 Other Inversion Techniques Birt les and More l (1979) propose an interpolation and inversion technique specifically suited for hydraulic head data that are sparsely scattered over a regional aquifer system. Also , the method can accommodate information such as borehole hydrographs, river and spring head measurements. The method involves superposition of a number of smaller problems. The method has elements of both direct and indirect inversion as indicated by Neuman (1973). The aquifer equation can be writ ten as: V • ( T V / i ) = S(^) - Q (1.35) which can be writ ten as V 2 / i ( x , y ) = -p(x,y) (1.36) where p(x ,y) is an unknown forcing function equal to — (l/T)\VT • Vh + Q — S(dhfdt)] The constraints on p(x, y) are such that values of h predicted at observation points should equal the measured values h*. Note that the above equation is linear, i.e.; the head (hi) at x, y can be constructed as: P = Pl+P2 + P3 + (1-37) and hi = h° + anpi + a t ) 2 P 2 + (1.38) where h® is the value of hi if p(x,y) = 0. Components p\, p 2 , etc., can be viewed as uniform values of infiltration received by sub-regions S\,S2, etc., and a^ y is the value of hi — hj. A matr ix system is thus formed: h = Ap + h° (1.39) 1.9 Other Inversion Techniques 22 and A is a N x M matr ix ( M is the number of mesh points). The matr ix A is filled with terms involving the difference between h{ from Poisson's equation due to a unit source at i and all other nodes. M a t r i x A is then reduced by deleting al l rows that do not correspond to observation points. This operation forms a new matrix Aq. A least squares solution for p can then be obtained as: which can be solved by standard algorithms. A set of interpolated values of h e s t at all points can then be found by multiplication of p wi th the original A matr ix: hest = Apest. Knowing hest at a l l points corresponds to the linear-direct problem, and the authors adopted a form of the inverse solution similar to Yeh et al . (1983) as described in the next paragraph. The authors concluded that after applying the method to an aquifer in England, computed patterns of T and S agree wi th the overall interpretation of the available data, and a new calibrated model based on the inverse values was capable of reproducing observed groundwater fluctuations to a high degree of accuracy. This method has not received much attention in the hydrogeologic literature. Yeh et al . (1983) used kriging to first obtain hydraulic head at all nodal points in the flow domain, and then used a direct method of inversion to obtain the inverse. In the algor i thm a biquadratic model of the drift in head was assumed and a kriged head map of the aquifer was made after calculating the variogram. A finite difference equation for the aquifer is wri t ten as: where A is the standard finite-difference global matr ix rewritten in terms of unknown transmissivity. This form of A is valid, provided that hydraulic heads are known at all nodal locations. In (1.41), n is the error estimated from kriging. In Yeh et al.'s paper, T is then estimated over P subregions by the Galerkin technique: p"« = [A^A^A^h* - h (1.40) A{h')T = b{h*) + r) (1.41) P M (1.42) 1.10 Overview of Chapters Two and Three G(i, j) is a basis function, and Te is the value of T estimated at the grid points. In matrix form: AGTe = b + v (1.43) and a solution is found by least squares: Test = (GTATAG)-lGTATb (1.44) T*st is the estimated vector of Te. The authors found that errors in the coefficients of the matr ix A depend on the kriged head values. They also demonstrated the use of the procedure on a synthetic example in which good results were obtained in spite of a high noise level in the data. For the problem studied, four values of T were the maximum that could be determined wi th the method. Alternate methods of inversion based on geostatistics and cokriging have been carried out (Ki tanidis and Vomvoris, 1983; Hoeksema and Ki tan id i s , 1985; Dagan, 1985). In the cokriging approach, a simple functional form for both the drift and the covariance of the parameters is adopted, thus yielding a small number of unknown parameters to be determined. The parameters are estimated through a form of a maximum likelihood or weighted least squares scheme. Then , a stochastic flow equation that relates hydraulic head and transmissivity perturbations is solved, followed by cokriging to yield the parameter estimates. Fundamental to this approach are the inherent assumptions regarding the form of the drift and variogram functions, and a data base in hydraulic heads and hydraulic conductivities (or transmissivities) sufficient to establish covariograms and variograms. However, in geotechnical problems (for example, stability of landslides, and design of open-pit mines) few hydraulic conductivity and hydraulic head data are usually available. As well , many of these problems are conceptualized on cross sections, wi th uncertain boundary conditions, and a drift (trend) in the hydraulic head data that cannot be resolved with a simple polynomial form. These facts make the application of geostatistical approaches difficult in problems of this k ind . 1.10 Overview of Chapters Two and Three 1.10 Overview of Chapters Two and Three A s mentioned in section 1.2, inverse goals can be categorized into three general forms. Given the type of inverse problem posed, a suitable inverse technique must be developed. The requirements of the technique should include: flexibility, robustness, and a minimum number of statistical assumptions either implici t ly or explicity implied in the method. It is the goal of Chapter 2 to introduce the idea of solving the inverse problem based on dif-ferent norms. A review of solution procedures is presented, and a construction formalism is adopted by introducing a series of objective functionals (norms) which are minimized to yie ld a variety of possible models. The inverse is carried out using a modification of a con-strained simplex procedure described by Box (1965) and Silva and Hohmann (1983). The algori thm requires no derivative computations and can be used to minimize an arbitrarily-complicated non-linear functional, subject to linear or non-linear inequality constraints. Linear constraints take the form of upper and lower bounds on the model, while non-linear constraints may take the form of (1.8). The algorithm produces a wide variety of acceptable models clustered about a global min imum, each of which generates data that matches observed values. Under the assumption that enough points have been sampled, the centroid of the sampled points is taken as the best estimate of the parameters. The inverse technique is demonstrated on a series of one and two-dimensional synthetic data sets, and on a sparse and inaccurate hydraulic head data set from Downie Slide, Bri t ish Co lumbia , Canada. Chapter 3 is devoted to developing the idea of a joint inversion scheme involving both thermal and hydrologic data. If data are inaccurate and sparse then model information may not be identifiable from the data. This situation ultimately effects the decision to adopt a particular inverse methodology. One way of overcoming data limitations (sparse hydraulic head measurements or few hydraulic conductivity estimates) in an inverse prob-lem would be to introduce an independently collected data set and apply simultaneous or joint inversion. The joint inversion method uses data from a number of different mea-1.10 Overview of Chapters Two and Three surements to improve the resolution of parameters which are in common to one or more functional relationship. One such data set is subsurface temperature, which is sensitive to variations in hydraulic conductivity (Smith and Chapman, 1983) in more permeable systems. The thermal and hydraulic head fields are linked by the fluid specific discharge. The best way to illustrate the advantages of a simultaneous inverse is through a simple example. Let us take the example of fitting a straight line through a group of data. Figure 1.1 is a graph of the equation y t = ax, + b. This equation (for i = 1 M , data points) can be wri t ten as: y = Fp (1.45) where p = (a, 6) and x2 1 F= . . (1.46) V XM 1 / and y is the theoretical predicted data. Let us now consider a non-unique problem where only one observed value of y* is available. A n infinite number of models (a, 6) can be constructed to fit a straight line through the point. The min imum length solution (see section 1.3, equations 8-10) in this case is: (0.5,0.5); however, this solution is not useful in predicting any other values of y{. Suppose that the model p was related to data by another function, for instance: Zj = a£j + b (j = 1,2,3, . . . .L) or z = Gp (Figure 1.2). If the observed data z* were perfect, the model p could be recovered exactly from this second expression, and knowing p, other values of yi can be predicted. The advantage of joint inversion lies in the usual case when data from both sets (y*,z*) are incomplete and inaccurate. Combining both data sets results in ( M + L) data wi th two unknowns for ( M + L — 2) degrees of freedom. A simultaneous inverse can be carried out by noting that: y* - Fp= e, z" -Gp = n (1.47) 1.10 Overview of Chapters T w o and Three These two equations can be combined into a larger-partitioned system: d* -Hp = $ (1.48) where (1.49) The least squares solution to (1.49) is: pest = ( H T H y l H T d * ( L 5 0 j Values of y at different locations, x, can be computed as: ypre = Fpest. A similar line of reasoning wi l l be adopted wi th hydraulic head and temperature data, wi th the exception that the two functional transformations involved are both non-linear, and the more powerful solution techniques presented in Chapter 2 wi l l be used. F i rs t , simple-hypothetical two dimensional examples of a joint inverse wi l l be carried out, followed by an example of a fictitious cross-sectional problem wi th an aquifer at depth. This example wi l l show the benefits of including thermal data if hydraulic head data are sparse. The inverse technique wi l l then be applied to the Downie Slide data; this time augmented wi th downhole thermal data. Chapter 1 - References References Backus, G . E . and J . F . Gi lber t , 1967. Numerical application a formalism for geophysical inverse problems, Geophys. J. Roy. Astron. Soc, 13, 247-276. Backus, G . E . and J . F . Gi lber t , 1968. The resolving power of gross earth data, Geophys. J. Roy. Astron. Soc, 16, 169-205. Bir t les , A . B . and E . H . More l , 1979. Calculat ion of aquifer parameters from sparse data, Water Resources Research, 15(4), 832-844. Box , M . J . , 1965. A new method of constrained optimization and a comparison wi th other methods, Computer Journal, 8, 42-52. Carrera , J . and S.P. Neuman, 1986a. Est imat ion of aquifer parameters under transient and steady-state conditions, 1, M a x i m u m likelihood method incorporating prior infor-mat ion, Water Resources Research, 22(2), 199-210. Carrera , J . and S.P. Neuman, 1986b. Est imat ion of aquifer parameters under transient and steady-state conditions, 2, Uniqueness, stabili ty and solution algorithms, Water Resources Research, 22(2), 211-227. Carrera , J . and S.P. Neuman, 1986c. Est imation of aquifer parameters under transient and steady-state conditions, 3, Appl ica t ion to synthetic and field data, Water Resources Research, 22(2), 228-242. Chang, S. and W . W - G . Yeh, 1976. A proposed algorithm for the solution of large-scale inversion problems in groundwater, Water Resources Research, 12(3), 365-374. Cl i f ton , P . M . and S.P. Neuman, 1982. Effects of kriging and inverse modeling on condi-t ional simulation of the A v r a Valley Aquifer in Southern Ar izona , Water Resources Research, 18(4), 1215-1234. Cooley, R . L . , 1977. A method of estimating parameters and assessing reliability of mod-els of steady state groundwater flow, 1. Theory and Numerical properties, Water Resources Research, 13(2), 318-324. Cooley, R . L . , 1979. A method of estimating parameters and assessing reliability of mod-els of steady state groundwater flow, 2. Appl ica t ion of statistical analysis, Water Resources Research, 15(3), 603-617. Cooley, R . L . , 1982. Incorporation of prior information on parameters into nonlinear regres-sion groundwater flow models, 1. Theory, Water Resources Research, 18(4), 965-976. Cooley, R . L . , 1983. Incorporation of prior information on parameters into nonlinear re-gression groundwater flow models, 2. Applicat ions, Water Resources Research, 19(3), 662-676. Dagan, G . , 1985. Stochastic modeling of groundwater flow by unconditional and condi-t ional probabilities: The inverse problem, Water Resources Research, 21(1), 65-72. Chapter 1 - References Delhomme, J .P . , 1979. Spatial variability and uncertainty in groundwater flow parameters: A geostatistical approach, Water Resources Research, 15(2) , 269-280. Dunbar , W . S . , 1979. Determination of fault models from geodetic data. P h . D . Thesis, Stanford University. Emsel lem, Y . and G . de Marsi ly , 1971. A n automatic solution for the inverse problem, Water Resources Research, 7 (5) , 1264-1283. Forsythe, G . E . , M . A . M a l c o l m and C . B . Moler , 1977. Computer Methods For Mathematical Computations. Prentice-Hall Inc., New Jersey. Freeze, R . A . , 1975. A stochastic-conceptual analysis of one-dimensional groundwater flow in non-uniform porous media, Water Resources Research, 11(45) , 725-741. F r ind , E . O . and G . F . Pinder , 1973. Galerkin solution of the inverse problem for aquifer transmissivity, Water Resources Research, 9 ( 5 ) , 1397-1410. Gambola t i , G . and G . Vo lp i , 1979. A conceptual deterministic analysis of the kriging technique in hydrology, Water Resources Research, 15(3) , 625-629. Glenn , J . E . and S .H . Ward , 1976. Statistical evaluation of electrical sounding methods, Par t 1, Experimental design, Geophysics, 4 1 ( 6 ) , 1207-1221. Glenn , J . E . , J . R y u , S .H. Ward, W . J . Peeples and R . J . Phi l l ips , 1973. The inversion of magnetic dipole sounding data, Geophysics, 38 (6 ) , 1109-1129. Hoeksema, R . J . and P . K . Ki tanid is , 1985. Comparison of gaussian conditional mean and kriging estimation in the geostatistical solution to the inverse problem, Water Re-sources Research, 21(6) , 825-836. Inman, J .R . , J . R y u , and S .H. Ward, 1973. Resistivity inversion, Geophysics, 38 (6 ) , 1088-1108. K i t an id i s , P . K . and E . G . Vomvoris, 1983. A geostatistical approach to the inverse prob-lem in groundwater modeling (steady state and one dimensional simulations), Water Resources Research, 19(3) , 677-690. Menke , W . , 1984. Geophysical Data Analysis: Discrete Inverse Theory. Academic Press, New York . Nelson, R . W . , 1968. In-place determination of permeability distributions for heterogeneous porous media through analysis of energy dissipation, Soc. Pet. Eng. J., 8 (1) , 32-42. Neuman, S.P., 1973. Cal ibrat ion of distributed parameter groundwater flow models as viewed as multiple-objective decision process under uncertainty, Water Resources Re-search, 9 ( 4 ) , 1006-1021. Neuman, S.P., 1980. A statistical approach to the inverse problem in aquifer hydrology, 3. Improved solution technique and added perspective, Water Resources Research, 1 6 ( 2 ) , 331-346. Chapter 1 - References Neuman, S.P., 1982. Statistical characterization of aquifer heterogeneities: A n overview. In Recent Trends in Hydrogeology, ( T . N . Narasimhan, ed.), G . S . A . special paper 189. Neuman, S.P. and S. Yakowitz , 1979. A statistical approach to the inverse problem in aquifer hydrology, 1. Theory, Water Resources Research, 15(4), 845-860. Neuman, S.P., G . E . Fogg and E . A Jackobson, 1980. A statistical approach to the inverse problem in aquifer hydology, 2. Case study, Water Resources Research, 16(1), 33-58. Oldenburg, D . W . , 1979. One dimensional inversion of natural source magnetotelluric ob-servations, Geophysics, 44(7), 1218-1244. Oldenburg, D . W . , 1984. A n introduction to linear inverse theory, IEEE Trans. Geosc. and Remote Sens., GE-22(6) , 665-674. Petr ick, W . R . , W . H . Pelton and S .H. Ward, 1977. Ridge regression inversion applied to crustal resistivity data from South Afr ica , Geophysics, 42(5), 995-1005. Si lva , J . B . C . and G . W . Hohmann, 1983. Nonlinear magnetic inversion using a random search method, Geophysics, 48(12), 1645-1658. Smi th , L . and D.S . Chapman, 1983. On the thermal effects of groundwater flow: 1. Regional scale systems, J. Geophys. Res., 88, 593-608. Tarantola, A . , and B . Valette, 1982. Generalized nonlinear inverse problems using the least square criterion, Ann. Rev. Geophys. Space Phys., 20(2), 219-232. Ward , S .H . , B . D . Smith , W . E . Glenn, L . Rijo and J .R . Inman, 1976. Statistical evalua-tion of electrical sounding methods, Part 2, Appl ied electromagnetic depth sounding, Geophysics, 41(6), 1222-1235. Wiggins, R . A . , 1969. Monte Carlo inversion of body-wave observations, J. Geophys. Res., 74, 3171-3181. Wiggins, R . A . , 1972. The general linear inverse problem: Implication of surface wave and free oscillations for earth structure, Rev. Geophys. and Space Phys., 10(1), 251-285. Yeh , W . W - G . , 1986. Review of parameter indentification procedures in groundwater hy-drology: The inverse problem, Water Resources Research, 22(2), 95-108. Y e h , W . W - G . , Y . S . Yoon , and K . S . Lee, 1983. Aquifer parameter identification with kriging and opt imum parameterization, Water Resources Research, 19(1), 225-233.3. Chapter 1 - Figures \ Figure 1.1 G r a p h of the line y = 0.3z+0.5 and min imum length solution to the inverse (y = 0.5x+0.5). Chapter 1 - Figures Figure 1.2 Graph of the curve z = 0.3f 2 + 0.5. C H A P T E R 2 Inversion of Hydraulic Head Data Summary A philosophy about the inverse is adopted that takes the form of constructing a variety of models under different norms. Only those features of the model that are consistently replicated are candidates for the ' true' model. In order to solve these problems, a power-ful opt imizat ion procedure is outlined which is based on the constrained simplex method. The algori thm easily handles constraints, does not require derivatives of the functional with respect to the parameters, and can be employed to minimize any arbitrarily com-plicated functional. The functional may involve more than one system of equations. The algorithm produces a cluster of points in parameter space about the global minimum of the functional. Each of these points generates data that matches observed values wi thin a prescribed tolerance. Under the assumption that enough points have been sampled, the centroid of the sampled points is taken as the best estimate of the parameters. The covariances of the estimated parameters and estimated data are easily computed. The inverse technique is applied to a sparse, inaccurate hydraulic-head data set observed at Downie Slide, Br i t i sh Columbia , Canada. In many inversion schemes, kriging techniques are employed to interpolate hydraulic heads to all nodal points in a finite-element scheme. It is found however, that on this cross-sectional problem with noisy data and a strong drift component, the interpolation is not possible. The difficulties are related to the inability of the residual kriging scheme to correctly identify drift in the presence of noise with a poorly distr ibuted data set. Four parameters are solved for: the free-surface position of the water table and three boundary conditions for the domain. Further simulations using 32 Summary a theoretical data set with assumed properties similar to that of Downie Slide show that wi th noise free data and an adequate spacing between measurement points, it is possible to interpolate an unbiased estimate of hydraulic head data at all nodes in the equivalent discretized domain. When the inverse technique is applied, the domain's hydraulic conduc-t iv i ty structure is correctly identified. The implications for Downie Slide are that in order to resolve more than a simple hydraulic conductivity structure on a cross-sectional domain, accurate field measurements of hydraulic head, a better spacing between measurements, and more detailed knowledge of the boundary conditions are required. 2.1 Introduction 2.1 Introduction The purpose of this chapter is to: • Introduce the idea of solving the inverse based on different norms. • Review some of the more popular optimization strategies used to solve the inverse. • Describe the inverse technique used in this thesis. • App ly the inverse approach at a field site. More specifically, the fundamental approach advocated in this thesis is to determine a vari-ety of different models of a hydrogeologic system, each of which generates data that satisfies the observed values. To do this, a series of different objective functions are minimized under a set of parameter constraints. A powerful optimization technique is detailed that can be utilized to minimize a complicated nonlinear functional without derivative computations. This inverse technique allows us to solve problems wi th more than one system of equations, or problems with uncertain boundary conditions, such as a free-surface. The method is similar to the Random Search Method (Price, 1977; Silva and Hohmann, 1983), and the Complex Technique (Box, 1965). The algorithm produces a simplex figure of an objective function (for example, equation 1.34) in parameter space by randomly distributing a set of search points wi th in constraint boundaries. For each search point, a value of the objective function is calculated. These search points are adjusted wi th in the constraints in a series of iterations unt i l the expected value of the objective function reaches a threshold value which is determined by the level of noise in the data. A statistical analysis involving all the sampled points is carried out when the parameter estimates are found. Wi th the technique described in this paper, the parameter space is randomly sampled to produce a covariance matrix; the only assumption being that the random sample be representative. The mean of ail samples is taken as the best parameter estimate. In addition, the covariance structure of the data is also provided, thus avoiding the effort of further conditional simulations. The algorithm only requires an init ial guess of the parameter values within a priori bounds. 2.2 Literature Review Theoretical aspects developed in this chapter are demonstrated with the inversion of synthetic data under different norms and an application of the technique to determining a model for the groundwater flow system in a landslide located along the Columbia River Valley in Bri t ish Columbia , Canada. The optimization strategy introduced here wi l l be applied to simultaneous fluid flow and heat transfer inversion problems in the next chapter. 2.2 Literature Review In the inverse approach measurements of hydraulic head, hydraulic conductivity or transmissivity, seepage flux, and the like are inputs to the inverse algori thm, and fitted hydraulic conductivity (or other parameters) become the output, along wi th the parameter covariance structure. Some authors (for instance, Clif ton and Neuman, 1982) use hydraulic conductivity values ( wi th their associated covariances ) determined from the inverse in conditional simulations. These authors showed that predicted hydraulic-head variances from the inverse/conditional simulations are greatly reduced over estimates of the variance obtained from kriging alone. Throughout this chapter the term 'model ' is defined as a configuration of material properties. Mathematical ly, the model takes the form of a vector m which consists of all unknown parameters (for example, hydraulic conductivities, boundary conditions) in a functional relationship that predicts physical data d (for example, hydraulic heads). B y a functional relationship, we mean some appropriate form of the groundwater flow equation. Current ly there are two broad classes of inverse methods documented in the hydrogeo-logic literature: a parametric approach, in which a continuous model is replaced by a finite number of parameters (e.g., Neuman and Yakowitz , 1979; Neuman et a l . , 1980; Neuman, 1980; Cl i f ton and Neuman, 1982; Carrera and Neuman, 1986a; Cooley, 1977, 1979, 1982, 1983; Yeh et a l . , 1983), and the cokriging-geostatisical approaches (e.g., Ki tanidis and Vomvoris, 1983; Hoeksema and Kitanidis , 1985; Dagan, 1985). Hybrids of both methods have also been developed (de Mars i ly et a l . , 1984). Fundamental to both of these broad 2.2 Literature Review cases of inverse techniques is the adoption of a model with a relatively small number of un-known parameters in relation to the number of available data points. The parameters are then determined through some form of maximum likelihood framework and/or a weighted least squares minimizat ion. In the cokriging approach, a simple functional form for both the drift and the covariance of the parameters is adopted, thus yielding a small number of unknown parameters to be determined. This approach requires a data base in hydraulic head and hydraulic conduc-tivities or transmissivities sufficient to establish covariograms and variograms. However, in many geotechnical problems (for example, stability of landslides, design of open-pit mines) it is usual that few hydraulic conductivity data are available and, as well , most of these problems are conceptualized on cross sections. Boundary conditions may also be uncer-tain. These facts make the application of geostatistical (cokriging) approaches difficult for these problems. Therefore, a geostatistically based method such as cokriging wil l not be pursued in this study. In the parametric approach it is typical to determine the parameters in a relatively small number of zones to ensure a mathematically unique solution is found, whether or not there is any physical justification for the parameterization (for example see Carrera and Neuman, 1985c, p 233-236; Birt les and More l , 1979). In some applications (for example, cross-sectional problems) the modeler may not possess sufficient geologic information to establish physically-based arguments for zoning sections of a finite difference or finite ele-ment grid. There are a number of possible methods for resolving this problem in particular cases. One approach is to allow the inverse method to identify zones of similar parameter structure. A n example of this approach wi l l be presented in section 2.12. In the solution of boundary value problems, one can face uncertainties in assigning a spatial representation for the hydraulic conductivity field, K , identifying sources/sinks and boundary conditions. This situation requires that an inverse scheme should be robust and capable of incorporating these types of conditions as unknowns. 2.3 Functional Relationships and Governing Equations 2.3 Functional Relationships and Governing Equations Solution of the inverse problem requires that both the forward and inverse problems must be clearly defined. A forward problem is set up by constructing a functional rela-tionship to predict physical data, given a set of input parameters to a physical model. The goal of inverse theory is to use a finite set of inaccurate observations to extract information about the model m (adapted from Oldenburg, 1984). For many physical problems the data and the model are related through a linear functional. However, in groundwater hydrol-ogy, when a parametric approach and a numerical scheme are used to solve the governing equations, a non-linear functional results: fe = 9f[i,m] (2.1) where 5 is a non-linear functional relating h, the data (the actual values of hydraulic head ), x, a vector of grid coordinates in a numerical scheme; and m , the actual model parameters, which could consist of hydraulic conductivities, boundary fluxes, sources and sinks. In applications of the inverse involving observed data, hydraulic heads have uncertain-ties associated wi th their values, resulting from interpolation or measurement errors. In these cases (2.1) takes the form: h* = S[x,m] + u (2.2) where / i ' is the data (interpolated or measured values of hydraulic head ), and u is a vector of residuals. Usually, v is assumed to have the following properties (a mean of zero and some covariance structure): E[y) = 0 E{uuT) = o\Vh (2.3) where a\ is a scaling factor used in defining the magnitude of the covariances. If kriging is used to determine V ^ , then a\ is taken as equal to one. 2.4 Construction In this study the functional S[x, m\ is given by the two-dimensional form of the steady-state groundwater flow equation: V K - V # = g (2.4) subject to - K - V * - n = 9/ (2.5a) on T i and * = g[x,y) (2.56) on T 2 . Here V = (d/dx,d/dy), V - = (id/dx + jd/dy), * is the hydraulic head, q represents sources and sinks, K is the hydraulic conductivity tensor, qj is a specified flux term on boundary Ti, and g{x,y) is a function specifying Dirichet boundary conditions on T 2 . When a numerical scheme is used to solve (2.4), the functional (2.1) takes the form: h = A-lf (2.6) where h is the approximate value of *i> due to the discretization, A is a global stiffness matr ix which is a linear function of K , and / represents a loading vector of the appropriate boundary conditions. 2.4 Construction Basic to most of the work in groundwater hydrology is the determination of a model which fits the data. The mathematical foundation for model identification has been devel-oped in a number of disciplines: control theory in electrical engineering (e.g. Schweppe, 1973), history matching in petroleum reservoir engineering (e.g. Gavals et al. , 1976), aquifer simulations in groundwater hydrology (e.g. Neuman and Yakowitz, 1977), and in geophysics (e.g. Parker, 1977; Oldenburg, 1984). A fundamental problem in determining a model which fits the data is that the 'true' model is a continuous function of the data. Unfortunately, there are an infinite number 2.4 Construction of functions that can reproduce a finite number of observations. In the present context, this could be called geologic non-uniqueness in that a decision must be made, based on the available geologic information, as to how to represent or parameterize the model. For instance, a decision must be made on the number of zones (layers) of similar hydraulic conductivity that exist in the flow domain. Work in this area has been carried out by Carrera and Neuman (1986a) using a max-imum likelihood criterion to aid in optimal selection of zones to which single values of the parameters are assigned. Generally, there is no unique method of resolving this type of non-uniqueness, unless strong geologic evidence exists to aid in the parameterization. As Oldenburg (1984) states: " This problem becomes more acute when data are inaccurate and sparse, and models which do not fit the data precisely are considered acceptable " . It is suggested that a reasonable approach to the inverse problem may be to construct a vari-ety of acceptable models which fit the data under different norms. These norms are based on the L\ and Lo lengths. This approach may provide insight into the non-uniqueness of the parameterization, and models which are geologically parsimonious, or the 'simplest' in some sense, can be selected as best approximations for the model. The underlying philosophy proposed in the chapter is to introduce a series of objective functions (norms) which can be used for model construction. The inverse can be stochas-tic or deterministic, depending on the assumed statistical nature of the interpretational model (see Appendix 2A) . In setting up an inverse problem, an interpretational model is formulated for which parameter estimates are sought. If insufficient geologic information exists to distinguish zones, a large number of parameters is incorporated in the model and the simplest geologic model can be chosen as the best estimate of the 'true' model. Simplicity can be defined as either the homogeneous-isotropic model, the smallest, the shortest deviatoric length from some estimate, the flattest, or the smoothest. The small-est model has the shortest Euclidean length of vector Y, and is the model obtained from (1.10) where the p^ term is set to zero. The S V D method described in Chapter 1 can be 2.4 Construction used to find p0. The smoothest and flattest models can be found by minimizing the first and second spatial derivatives, respectively between physically adjacent model parameters (see Menke, 1984). Assuming that the data are free of gross errors, a difference in model parameters estimated using these different norms suggests a possible inconsistency in the model structure. Conversely, differences in constructed models may also point to gross errors in the data if the parameterization is correct. Silva and Hohmann (1983) used such an approach in treating a non-linear magnetic inversion problem in geophysics. The idea of using different norms for construction purposes is well established (see Oldenburg, 1984 for a review). T w o simple examples wi l l be presented in the next section to illustrate this approach. Generalized Norms A generalized Li norm is introduced as (after Neuman and Yakowitz, 1979): J = (h - h~)TV-l{h - h") + X(Y - V)TV-l(Y - Y') (2.7) where the above terms are defined as (Neuman and Yakowitz , 1979): and Vy are head and log conductivity covariance matrices, Y is a vector of log conductivities determined by the inverse method, Y* is a vector of observed or estimated log-hydraulic conductivities, h" is a vector of observed or estimated hydraulic heads, and A is a scaling factor, which may be unknown. and Vy are defined based on the characteristics of the data set and the finite element mesh. If h/ and Y' are estimated by kriging, along with Vh and V y , then A = 1. Wi th an unknown A we recognize that we may have knowledge about the structure of the covariance matrices but not their magnitudes. Equat ion (2.7) can be derived from a maximum likelihood consideration for a Gaussian distribution (Carrera and Neuman, 1986a), but one does not have to assume any underlying statistical distributions of h or Y to apply the norm (Schweppe, 1973). The objective function can also be viewed as a weighted sum of Li prediction error (heads) and L 2 solution simplicity (log K ) . The optimization technique outlined later in this chapter is capable of minimizing 2.4 Construction complex functionals, and although not all the listed functionals are used, some forms of (2.7) are presented here for completeness. If no prior information is available on log-conduct ivi ty the following norms can be used: 'Smallest model ' J = (h- h'fV-'ih-h*) + A || Y \\l (2.8) Subject to: Yi < Y < Yu where || • 11^ = 5Zj = i ( ^ j ' ) 2 - ^ ' a n ^ Yu are prior lower and upper bounds on log-conductivity. Equat ion (2.8), when linearized (for example in a quadratic programming scheme) reduces to the more familiar ridge-regression technique (Draper and Smith, 1981). Note that under this scheme, lower log-conductivity values are penalized, such that minimization of the norm finds the highest log-conductivities (within the constraints) that satisfy the data. As noted by Menke (1984) this norm is an inappropriate measure of structural simplicity, in that it tends to enhance variabili ty between parameter estimates assigned to different zones. Other forms of this norm, such as the flattest or smoothest are preferred for this purpose. One could also minimize the length away from some assigned value: 'Smallest deviatoric model ' J = {h- h')TV-l[h - h~) + A || Y - y \\l (2.9) Features that exist in these models are smoothed artifacts of the real system. Another group of objective functions under the L\ norm are suitable for construction purposes. This norm is generally considered to be more robust than the Li norm as the underlying distr ibution in u is viewed as being drawn from an exponential probability 2.5 Example Construction Problem distr ibution function (pdf). These pdf's have a larger portion of their tails extending outwards from the centroid than does a Gaussian distribution. Consequently, outliers in a data set, due to systematic and/or gross errors, or incorrect parameterizations, do not effect the norm calculation to as large a degree as do Li norms (Menke,1984). The drawback in applying the L\ norm is that statistics of the parameter estimates are not readily available. Some of these objective functions are: J = \h - h'\ + X\Y - Y"\ (2.10) or a weighted L\ norm: J = yU il + xy1-^ «J (2.11) ]=1 ' 1=1 where o^. and oyt are the standard deviations of an observation of h'^ and Yt~, respectively. 2.5 Example Construction Problem Figure 2.1 shows an example of an inverse problem which involves fitting a straight line through a set of data points. The data points are generated using the quadratic equation y = 1.0 + 0 . 5 i 2 . The ' true' model in this case consists of two parameters m = (1.0,0.5). Four noise free values of y form the data base. If it is assumed that the underlying interpretational model is a linear equation, a unique least-squares fit [Li norm) yields an inverse of m = (0.31,1.96), while the L\ solution is m = (1.0,1.47). Although these two results are unique mathematically, the inconsistency in the results points to an incorrect interpretational model. In this case the modeler may be rightfully concerned that the underlying 'physics' of the process may not be correct. A second example is shown on Figure 2.2 which shows a one dimensional steady-state groundwater flow system. The 'true' model consists of four layers of equal thickness but w i t h different hydraulic conductivity values. The boundary conditions consist of a specified head at one end (*(/) = hi), and a specified dimensionless flux (*'(0) = q*) at the other. 2.5 Example Construction Problem 43 The hydraulic head solution to this problem is: i rff + q'K{0) i (2.12) Note that for the boundary conditions that are imposed, and having 9 values of hydraulic head at the locations shown on the figure, the inverse problem is non-linear. The problem can be transformed into a linear one by recognizing that only one set of ratios in hydraulic conductivity for the various layers wi l l produce the observed hydraulic head pattern. Let us assume that the correct interpretational model is made, in that we know the location and thickness of the four layers present. The layer thicknesses are constant, and are set equal to 2dz. Mul t ip ly ing the integral expressions in (2.12) by K(0) yields: where Px = j^}, Pi = 7^(3)' a n ( ^ ^ 3 = 7c(4) Note that K\ — A'(0) in this problem. Now, the four layer hydraulic conductivity values are reduced to three ratios [P\, Pi, Pz)- The ratio in the first layer is fixed at 1. If (2.13) is written at the nine observation points as shown in Figure 2.2, a set of simultaneous linear equations results: where m = (P\, P2, P%) and G is a matr ix containing the geometry of the problem. A series of three simulations is presented. Hydraulic conductivity values assigned to each layer are 10°, 1 0 _ 1 , 1 0 ° , and 10~ 3 m/sec, respectively. Hydraul ic head values are obtained for each of the nine data points by first solving the forward problem. The boundary conditions are: hi — 1.0 and q" — -1 .0 . In one inverse simulation, using the 9 noise free values of hydraulic head, identical solutions are obtained using the L\ and 1/2 norms. In the second simulation a gross error is introduced into the eighth data point [hg). Table 2.2 shows the inverse results based on this data; note the degradation of the model under the L2 as opposed to the L\ norm. Note that the L\ norm algorithm (2.13) d = Gm (2.14) 2.6 Optimization Technique fit 2 data points exactly and ignores the remaining data points. In the third simulation an incorrect interpretational model is chosen. In this case a two layer model is chosen to represent, in reality, the four layer system. A n inverse based on the original 9 values of hydraulic head (Table 2.1) shows an inconsistency between model parameters (Table 2.3). Therefore, assuming the data are blunder-free, a difference in two constructed models based on different norms may indicate an incorrect parmeterization. Geologically based arguments would then have to be used to reevaluate the parameterization. 2.6 Optimization Technique Equation (2.7) represents a^constrained non-linear optimization problem and a large number of optimization techniques can be used to solve it. Some algorithms use quadratic programming techniques to solve a linearized form of (2.7). These programming techniques involve the solution of a set of simultaneous linear equations, and are sensitive to rank deficiency. Rank deficiency occurs when one or more of the unknown parameters are linearly dependent. This situation can occur if a large number of unknown conductivities are estimated in an overparameterized model with inadequate prior information. If rank problems are encountered in a conventional quadratic programming scheme, the matrix solution routines fail to yield a solution. Cooley and Naff (1986) describe instances of this type of problem. Other optimization techniques involve the minimizat ion of an objective function in its original , non-linear form. Some methods require that derivatives of the objective function wi th respect to model parameters be calculated, others do not. The derivative or gradient based methods as they are called are iterative in nature. These methods are usually faster than non-gradient methods but may diverge (fail to find an answer), or converge to a local rather than a global min imum (Reklait is , et a l . , 1983). A non-linear inverse problem may have a complex functional surface, and there is no guarantee that any technique will converge to a global min imum. Gradient methods work best if the inital guess is linearly 2.6 Optimization Technique close to the solution. To investigate non-uniqueness it would be desirable to repeat the procedure many times with different init ial guesses. Unfortunately, with gradient based methods it is sometimes difficult to arrive with a new initial guess for the parameters, especially if the problem involves constraints. Carrera and Neuman (1986b) evaluate sev-eral gradient based methods and show examples of non-convergence using examples in groundwater hydrology. M o d e l covariances are useful end-products of an inverse scheme. However, because the dis t r ibut ion for mest may be non-Gaussian, the covariance of the estimated model parameters may be difficult to interpret, especially in terms of confidence intervals. The covariance of the parameter values can be estimated as (Menke, 1984): [cov mest\ = G~9\cov h'\G-gT + [I - Rn]\cov m'}[I - Rn\T (2.15) where Rn — G~9Gn, I is the identity matrix, and the final n'th iteration is used. Also: [Gnliy - a m j (2-16) G~g is called the generalized inverse of G (see Chapter 1 for details). The reason for the non-Gaussian nature of the J surface has been examined at a theoretical level by Tarantola and Valette (1982). Because (2.7) is non-linear the probability distribution function of J w i l l in general be non-Gaussian, although fortunately the experience in groundwater hydrology wi th aquifer inverses is of Gaussian or near-Gaussian behaviour (eg. Cooley, 1977). Nevertheless, the maximum likelihood point of a non-Gaussian objective function may not yield the most sensible parameter estimates. Gaussian distributions are symmetric , so the max imum likelihood point always coincides wi th the mean value. For an arbitrary non-linear surface (for example, multi-modal, skewed ) the maximum likelihood point can be quite far from the mean value. Figure 2.3 shows schematic diagrams of objective functions for several non-linear problems. Figure 2.3a (after Menke, 1984) shows an objective function that has several local minima, or points of non-convexity with a well defined global min imum. Plots (b) and (c) of this figure show problems with two or more 2.6 Optimization Technique solutions, and (d) shows a fiat-bottomed objective function wi th a finite range of solutions. Plot (e) shows a problem with an infinite range of solutions. It is the goal of inverse theory to condit ion the norm in such away that wil l yield features like (a), wi th a well defined solution. The technique proposed in this chapter samples the objective function surface close to the minimum to find an average value of the model parameters. These values may be more appropriate than a single estimate obtained by most algorithms. The inverse procedure produces estimates of the data as: hest = %(x,mest) (2.17) where hest is the estimated hydraulic head data produced from the inverse-parameter estimates. Ordinarily the correlation structure of hest cannot be obtained from the inverse approach even though the prior covariance structure of the data may be known. In order to generate this matrix (cov hest) one has to use the Monte Carlo approach (Clifton and Neuman, 1982), or use a perturbation technique (Townley and Wilson, 1983) after obtaining the estimated model covariances, (cov mest). In order to deal with solution divergence and false convergence, and to avoid having to compute gradients, inflexibility of starting guesses, and linear approximations in covariance calculations, we can adopt alternate methods of solution. Silva and Hohmann(l983) used a controlled random search technique to optimize a non-linear functional in magnetic anomaly interpretation. In this thesis we use a procedure that is found to be similar but more efficient than theirs. A constrained simplex approach is used (Complex technique: Box , 1965) for global optimizat ion. The algorithm easily handles constraints, does not require part ial derivatives of the functional with respect to the parameters, and can be employed to minimize complicated non-linear functionals which may contain more than one system of equations. In addit ion, convergence to a global minimum is enhanced (even in the presence of local min ima ) provided the global minimum exists within the parameter constraints. Besides finding the global minimum, the algorithm produces a sampling of the parameter space around the global minimum. This sampling produces the parameter and 2.6 Optimization Technique estimated data covariances matrices. In cases like Figure 2. 3d upper and lower bounds to the model may be resolved wi th this method. The constrained simplex technique is well known in the chemical engineering and op-t imization literature (Reklaitis et a l . , 1983). The algorithm consists of choosing L search points (L > N + 1 parameters) randomly within the constraint boundaries. This initial task can be accomplished wi th a random number generator. These points set up the initial simplex figure. Each search 'point ' is a vector of model parameters. Thus , L points form an N x L matrix P". The objective function is then evaluated at each search point and the search points with the min imum and maximum objective function values are identified. The parameter estimates and function value for each point are stored in computer memory. A t the start of each iteration, the point with the maximum function value is reflected and expanded about the centroid of the remaining points according to the relation: B = S + a(S-Q) (2.18) where B is a vector of new parameter values, S is a vector of parameter values with each entry calculated as the average of the parameter estimates using the other (L — l ) search points, Q is the previous position vector of the maximum J value, and a is an expansion factor. Box (1965) recommends a value of a = 1.3, based on experience. If the reflected point violates a constraint for one or more of the parameter values it is relocated just slightly within the appropriate constraint boundaries. The objective function is then reevaluated at B. If the point B persists in having the worst function value it is contracted one-half way towards the center of the remaining points. If the point B has a smaller objective function value than before, then new minimum and maximum values are sought. This series of steps constitutes one iteration. In this way the centroid of all the values proceeds to the global minimum, and is not affected by local min ima or saddle shaped troughs in the objective function surface (Schwefel, 1981). If the min imum of all tr ial points remains the same for a number of iterations (usually set to 100), then the entire simplex figure is moved one half way towards the min imum. This modification to the 2.7 Stopping Criteria constrained simplex procedure is made to speed convergence. The iteration procedure is terminated when all search points have objective function values less than a specified noise level. Under ideal (noise-free) conditions the objective function (2.7) behaves as a quadratic surface and the min imum objective function value corresponds to the true parameter values. However, when the data are contaminated by noise, the min imum does not correspond to the true parameter values. Its expected value is equal to the noise level in the data. This behaviour is discussed in more detail in the next section. 2.7 Stopping Criteria In the presence of noise, it is important that the data predicted by the inverse does not reproduce the observed data exactly. A model generating a data misfit below its expected value wi l l show structures that are artifacts of the data. The expected value of the data misfit is related to the noise in the data set. This point is discussed in the next paragraph. In the converse situation, if the data misfit is too large then information about the model contained in the data wi l l be lost (Schlax, 1984). In the constrained simplex procedure, the centroid of the L search points moves towards the global min imum. W i t h each succes-sive iteration the entire simplex figure becomes more compact. Ult imately, the centroid would move to the min imum and the entire simplex figure would collapse to the centroid value at that point. As mentioned, this maximum likelihood point may not give the most sensible estimate of the parameters in all cases. As well, meaningful covariances could not be generated. Conversely, it is undesirable to underfit the data and arrive at poor covari-ance estimates. Therefore, an appropriate stopping cri teria must be defined. As noted by Fullagar and Oldenburg (1984), the appropriate level of misfit is a subjective choice. Note that no assumptions are made regarding the underlying statistical distribution of the functional, the parameters, or the data when the parameter space is sampled to compute covariances. However, in order to terminate the optimizat ion procedure, it is necessary that some assumptions be made in this regard. When working wi th Li norm functionals, 2.7 Stopping Criteria 49 the assumption is made that the values of the objective function are normally distributed about some mean value. In the constrained simplex algorithm, when the average value of the L objective function values is below the expected value of noise, the procedure is terminated and a sampling of the parameter surface is carried out using all search points. The appropriate measure of misfit is the x 2 test when data errors are Gaussian, inde-pendent wi th zero-mean. Here the x 2 misfit is defined as: X 2 = D * (2-19) i=i and di is a computed data point. In the controlled random search method / is replaced by M. If J > 5 the expected value of x 2 is approximately / and the standard error of x 2 is oxi = y/21. Any model is defined as acceptable when / - axi < x2 < I + o-x? (2.20) The most likely model is the one corresponding to x 2 = I- Note that objective function (2.7) can be writ ten in the form of (2.19) through the use of a spectral (eigenvalue) decom-posit ion. Therefore, we may apply this statistical test to (2.7). Tn (2.7) it is recognized that the expected value of the functional (noise level) is equal to o\M, because the median of X 2 is approximately equal to the number of degrees of freedom, and the objective function has been impl ic i t ly scaled by a factor of o\. If the L\ norm functionals are minimized ( equations 2.10 &: 2.11) then a different strategy is adopted. These functionals can be viewed as following exponential distributions. For the L\ norm the expected value of x 1 is ( j f ) a M wi th a variance of (1 - \)M (Parker and M c N u t t , 1980). 2.8 Computat ion of Covariances 2.8 Computation of Covariances The constrained simplex algorithm produces L acceptable models, each of which satifies the data, and the centroid of the global minimum is expected to be more meaningful then the single estimate obtained by most algorithms when the surface of the objective function is non-Gaussian. This approach produces equivalent results to a Monte-Carlo solution to the inverse (see Parker, 1977). Therefore, the average of all the L search points is taken as the best estimate of the parameters. The parameter space can also be sampled to obtain the parameter covariances, without assuming linear behaviour close to the min imum value. An advantage to this strategy is that covariances of model parameters may be established regardless of the norm that is minimized. However, assumptions on the statistical nature of the covariances must be made in order to compute confidence intervals (for instance, Gaussian for L 2 , exponential for L i ) , and to choose a stopping cri terion. Silva and Hohmann (1983) adopt this strategy. Wiggins (1972) also showed that model covariances can be recovered from repetitive Monte-Carlo samples of parameters. These covariances are estimated: c o v m ? = JT=T) D ^ * ~ Pil\pjk - Pj] (2-21) where cov m^1 is an element of the covariance matrix, P*k is the k ' th estimate of the i 'th parameter and, p\ is the average value of the i ' th parameter such that: 1 L Pi = I E pik (2-22) k=l A n advantage of the method is that we may also compute the covariance of hydraulic head estimates at data points, (cov heat) by storing the computed hydraulic head values in a mat r ix H' that correspond to each parameter estimate. For each search point (k = 1....L) there wi l l be a vector of length M of hydraulic head values. Then, 1 L [cov he%3 = ^—T) - h*}[H;k - (2.23) 2.9 Selection of Xopt 51 and L * i = T (2.24) ^=1 In the results section, only the standard deviations in hest w i l l be examined. We do not consider interpretation of the off-diagonal terms. 2.9 Selection of Xopt For an extensive discussion of the determination of Xopt the reader is referred to Neu-man and Yakowitz (1979). Similar rules can be adapted for our particular optimization scheme. 1. o\ and o\ are known. Then a2 ^oPt — ~y (2.25) In the case where [cov h"\ and [cov m'\ are estimated from the data, by say kriging, then one would set A = 1. 2. Either o\ o r . a y is unknown or both are unknown. In the case of an assumed [cov m*] as in the Bayesian type or a Fisher type (see Appendix A) where little is known about the absolute level of noise in the data, one can proceed wi th a 'Stagewise Opt imiza t ion ' approach (Carrera and Neuman, 1986a). i) Select A n = 1 for n = 1 and solve as in 1 above. ii) Compute 2 (Y — Y~)TVyl(Y — Y") o\ = K- >—JLA L (2.26) 2 _ {h-h')TV~\h-h^) M where o\ and o\ are estimates of the variances in log-conductivity and hydraulic head. Compute * — f t 2.10 Summary of Algori thm iii) Repeat steps (i) and (ii) unti l convergence in A is noted. As noted by Carrera and Neuman (1986a), either c\ or o\ can be treated as unknowns and solved for in the optimization procedure. In some cases a functional like (2.7) or (2.11) can be minimized for a succession of A values and the most appropriate value of A chosen graphically. 2.10 Summary of Algorithm In this section the computational strategy used in solving an inverse problem wi l l be summarized. 1. Decide on an appropriate interpretational model and parameterize the flow domain. Assign boundary conditions that are assumed to be known or treated as parameters to be solved for in the inverse. Discretize domain for the numerical solution to the boundary value problem. 2. Choose the functional (norm) to minimize, for example: J = (h — h*)TV-l(h - h*) + X{Y - Y')TV-l[Y - Y*) (2.27) This can be written as: J = Jh + XJY (2.28) w i th constraints: Yi<Y< Yu where Yi and Yu are lower and upper bounds on the log-conductivity parameter esti-mates. 3. Set the random number seed and randomly pick L t r ia l points wi thin the constraint boundaries. Each trial point is a vector of N parameter values. 4. Compute the hydraulic head at each nodal point for each of the L t r ial points. 2.10 Summary of Algorithm 5. Compute the functional J for each set of trial Y values. The functional J f t can be computed by Cholesky Decomposition as: Jh = (h- h~)TV-\h - h') = vTV~1v (2.29) Let V = V^1v (2.30) or Vhy = v (2.31) so that if a L L T decomposition of Vk is used: L L T y = v (2.32) where L is a lower triangular matrix. This step is followed by two back substitutions to obtain y, then Jh = vTy (2.33) A similar procedure is followed for Jy. Note that the inital decomposition is carried out only once. For each functional evaluation, two back substitutions and a dot product are required. 6. Follow the constrained simplex procedure until the global minimum is reached at the appropriate noise level. 7. Compute the covariance of h e s t and mest. 8. Check for non-uniqueness by resetting the random number seed and resolving the problem. In concluding this section a summary is presented of the major advantages and disad-vantages of the method: Advantages 1. F l ex ib i l i t y in the form of the functional. 2.11 Inversion of Theoretical Data 54 2. Covariance of the data cov hest is computed as part of the inverse solution. 3. Less sensitive to computational instability. If inverse fails (i.e. global min imum does not exist), repeat the procedure for a different norm. 4. Easy to generate different initial points to establish uniqueness of solution. 5. Programming ease. Disadvantages 1. A n extensive number of repetitive computations must be performed. D a t a on compu-tational aspects wi l l be presented in the results section. 2. Because of computational demands, practical application is limited to steady state problems at this time. 2.11 Inversion of Theoretical Data In order to investigate the robustness of the proposed approach, the technique is applied to a succession of one and two-dimensional problems. In some cases Gaussian random noise is added to theoretical data to simulate corruption of 'true' data values. Alternatively, errors in the data could contain structure, typical of kriging estimates. Only the former case is analyzed in the theoretical examples presented here. Various statistics are computed wi th each simulation. Two of these statistics yield an estimate of the,errors in the data. These statistics are: y=i The following two statistics give a quantitative measure of how well the inverse reproduced the actual values of the model and data (Yact and hact): (2.34) (2.35) (2.36) j=i 2.11 Inversion of Theoretical Data i M S& = ( ^ D ^ - f c i C T ) * (2-37) j = \ where Y and h are the centroid of sample values of hydraulic head and log-conductivity from the inverse algorithm, respectively. One Dimensional Examples In the first series of synthetic examples, we examine the ability of the inverse method to reconstruct a realization of a log-conductivity field in a one-dimensional, discretized domain . Firs t , a sequence of log-conductivity values is generated along a one dimensional gr id of unit length, using a multivariate random number procedure (e.g. Clif ton and Neuman , 1982, p 1219). The grid is regular and consists of 19 elements and 20 nodes. A n exponentially-decaying covariance function is adopted to describe the spatial continuity of the hydraulic conductivity variations. Various log-conductivity data sets are generated from a distribution having normalized (with respect to the length of the grid) integral scales of (-7 ) 0.8 and 0.2, and standard deviations (ay) of 0.1 and 1.0. The integral scale measures the distance at which the correlation decays to a value of e _ 1 . The mean value of log-conductivity E(Y) in all cases is -2.0. The various realizations are shown on Figure 2.4, and are grouped into categories A , B , or C as shown on Table 2.4. A one dimensional form of the finite element system (2.6) is solved wi th boundary conditions h=0.0 at one end and h=1.0 at the other. For each realization the actual hydraulic heads at each node and actual Y values within each element are corrupted with Gaussian random noise at various percentages of the true value. In all cases Vh and Vy are taken as diagonal matrices with known terms assigned from the noise corruption procedure. In this series, no zones or groups of elements are chosen; the inverse is free to assign hydraulic conductivity values in each of the finite elements. Results are given in Table 2.4. For one simulation in each of the three groups (Table 2.4), a priori hydraulic conductivity information has not been incorporated in the inverse. Note that the inverse under these circumstances is non-unique. Upper and lower bounds on the Y parameters are set to -4.00 2.11 Inversion of Theoretical Data and -0.1, respectively. For example, (Table 2.4, line 1) the first simulation uses hydraulic head data that is corrupted wi th Gaussian random noise having a standard deviation equal to 10 % of each value. This corruption leads to an S^- noise value of 0.074. Using parameter estimates from the inverse, the actual hydraulic head values are reproduced to within a value for of 0.014, clearly less than the noise level of 0.074. The results for the L2 norm (equation 2.7) without the penalty term) show in general that the inverse procedure reduces hydraulic head variances over the noise in the data. The actual log-conductivities are reproduced to a Sy value of 0.424. However, for this cases the actual values of Y are not replicated that well (lines 1, 4, and 7). This behaviour is considered to be due to the non-uniqueness inherent in the type of boundary value problem being solved. For the remainder of the runs in each grouping penalty term in (2.7) is applied, incorporating the noisy a pr ior i Y' information. For example, incorporating prior hydraulic conductivity data in group A with a noise level of 0.20 leads to a standard deviation in the estimated values of hydraulic conductivity of 0.152; or a reduction of 25% in comparison to the noise level (line 2). In the other trials the reduction varied from about 1% to 42%. Even the more complicated realizations (<7y= 1.0 : Group A) are resolved to suprisingly good detail . Two Dimensional Example In order to further verify our solution approach, a comparision is made wi th a solution recently published by Carrera and Neuman (1986c). The example is a hypothetical square aquifer of 36 k m 2 , subdivided into 9 zones wi th transmissivity values as shown on Figure 2.5. Eighteen observation points are located within and along the boundaries of the aquifer. Boundary conditions are shown on the figure and are discussed in Carrera and Neuman's paper. In one series of simulations these authors first solved the steady-state hydrologic problem, and then corrupted the hydraulic head data wi th Gaussian noise of various levels. Next , they applied their inverse procedure both with and without prior information. Our solution obtained with the constrained simplex technique is compared to the solution 2.12 Analysis of Field Data presented by Carrera and Neuman for the steady-state situation without prior information. The steady-state problem is solved by subdividing the problem area into 36 square elements (72 triangles) and applying the finite element technique (equation 2.6). The inverse problem consists of determining one unknown log-transmissivity value in each of the zones, using the same zoning pattern as Carrera and Neuman (1986c). Each zone contains 8 triangular elements. The actual hydraulic heads at the 18 measurement points are corrupted by Gaussian random noise prior to solving the inverse problem. Two realizations are chosen wi th different random seed numbers. B o t h trials yield realizations of N(0,1 m). The boundary conditions are assumed to be known wi th certainty and the L2 norm functional (without the conductivity penalty term) is adopted to solve the inverse problem. The lower and upper bounds on the log-transmissivity parameters are set to 0.1 and 2.5. The number of search points is 30. Typ ica l runs involved 5000 iterations and took about 10 minutes of C P U time on a F P S 1 6 4 / M A X array processor. The results of two of our runs and one of Carrera and Neuman's are shown in Table 2.5. The sum of the squared errors (2.36) for two of our trials and three of Carrera and Neuman's are also computed. Because of the random populations from which the noise is drawn, we are not comparing exactly the same data as did the previous authors. However, the results suggest that our procedure reproduces the actual transmissivity field at least as accurately as Carrera and Neuman. 2.12 Analysis of Field Data In this section we discuss the results of applying the inverse theory presented in this thesis to a field problem. The site chosen is the Downie Slide in the east-central part of the Province of Br i t i sh Columbia , Canada. The slide has been studied in the period from 1965 to the present time. Relevant data on Downie Slide has been presented by Piteau et a l . (1978), Imrie and Bourne (1981) and Patton (1983, 1984). Previous groundwater modeling of the Slide was carried out by Freeze (1977, 1982) and Dalton (1986). 2.12 Analysis of Field Data 1) Hydrogeology of Downie Slide Downie Slide is a 1.5 x 10 9 m 3 landslide located on the right bank of the Columbia River about 60 k m north of the Revelstoke D a m . The slide has been studied extensively since 1973, and from 1978 to 1982 a major drainage program was carried out which included installat ion of 2,432 m of drainage adit and 13,720 m of drainpipe to offset increases in pore-pressures caused by the flooding of the toe of the slide by the Revelstoke Reservoir. The slide is believed not to have moved quickly or catastrophically, and probably moved in large intact blocks. Most of the movement is believed to have occurred on a major shear plane which runs at about 250 m depth below surface. The basic stratigraphy of the slide is well known on the basis of surface geologic map-ping, adit mapping and borehole logging. The rock formations in which the slide is situated are composed mainly of schists and gneisses of the Shuswap Metamorphic Complex. The rocks are highly foliated, and with the exception of projecting some of the larger shear zones, hole-to-hole geologic correlation is not possible even between closely-spaced bore-holes (20 m ). Sources of hydrogeologic data consist of the following: regional and local geologic mapping, geological and geophysical logs obtained from boreholes, a large-scale recovery test carried out in one of the adits, piezometer response testing, hydrological and meteorological data gathered at the site, discharge records obtained from two adits, and computer simulations of surface water and groundwater hydrology. In order to set geometries and boundary conditions for the inverse problem at this site, an examination of regional hydrogeology is necessary. This examination is carried out by modeling a geologic cross-section which represents the regional flow patterns from major axes of symmetry located at the Revelstoke Reservoir, and the mid-point of the Monashee M o u n t a i n range (Figure 2.6). The section is modeled as a homogeneous, isotropic system wi th a specified fluid flux across the water table. This surface is expected to be close to ground surface. Recharge along the section is provided mainly by snow melt in the summer months. F low in the section is approximately horizontal due to the length/depth ratio of 2.12 Analysis of F ie ld Data the flow system. Most of the discharge occurs through the adit system (modeled as a point sink of constant head), with a minor amount of flow to the Revelstoke Reservoir. In this study, we focus on the section bounded by C — C on the west and the Revelstoke Reservoir on the east. Figure 2.7a shows a plan view of the Slide wi th drillholes completed since 1965. Also shown are some of the principal intrumentation locations. Figure 2.7b shows a cross sectional view of the slide. None of the aforementioned tests, studies, or drainage experience have indicated any coherent stratigraphy on the slide. There is no evidence to suggest that the hydraulic conductivi ty varies to a large degree either above or below the major shear plane. The shear zone itself probably behaves locally as a low hydraulic conductivity barrier, however on the scale of the slide and the mathematical solution this feature can be neglected. It is expected that groundwater movement wi l l be controlled primarily by fracture pathways. Given the scale of the analysis, the rock mass is modeled as a continuum. Freeze and Cherry (1979) report hydraulic conductivity values of 10~ 8 to 10~ 4 m/sec for this type of mater ial . Results of the various hydraulic conductivity estimates are summarized in the next section. 2) Hydraulic Conductivity Distribution Estimates of hydraulic conductivity determined at the slide come from four main sources: borehole measurements during dri l l ing, piezometer response tests, a large-scale recovery test, and computer-model calibration. A l l of the values obtained from borehole tests dur ing dri l l ing have been discounted due to heavy layering of bore walls wi th dr i l l mud and other additives. Response tests carried out in piezometers indicate a range of hydraulic conductivity values from 1 0 - 9 to 10~ 6 m/sec. However, in reviewing the reli-abili ty of these tests, Freeze(l982) discounted many of the individual results because the tests were carried out over short terms and in long portions of the boreholes. T w o of the more reliable sources of hydraulic conductivity data come from the recovery 2.12 Analysis of Field Data test and computer-model simulations. In the recovery test, two flowing adit holes were closed off and the response was measured in two piezometers at large distances away from the adit. Meneley (1977) computed a value of 6.2 x l 0 ~ 6 m/sec for the hydraulic conductivity of the slide mass above the shear zone. Freeze (1977) carried out a suite of steady-state simulations of groundwater conditions on the slide. Cal ibrat ion runs showed that on the basis of the observations at that time, a homogeneous-isotropic flow system with a hydraulic conductivity of 10~ 7 m/sec could not be distinguished from the flow system for a layered medium (10~ 6 m/sec for the rock mass and shear-zone hydraulic conductivity of 1 0 - 8 m/sec). Addit ional simulations and drainage experience provide further support for a homogeneous, isotropic configuration (Freeze, 1982). 3) Recharge Rates and Water Budget Recharge in the slide area is considered to be concentrated in a 3.2 x l O 6 m 2 area at the top end of the slide. The evidence to support this conclusion comes from two main considerations: 1. Surficial deposits in the area of the head scarp consist of blocky rockfall of high porosity. This material would allow for almost complete infiltration. 2. The surficial deposits on the lower two-thirds of the slide consist of a th in veneer (up to 10 m thick) of glacial t i l l that would tend to promote surface runoff. Recharge to the slide is estimated to be 9.1 x 10 5 m 3 / y r (Freeze, 1982). Over the recharge area, a flux of 3.0 x 1 0 _ 1 m / y r can be calculated. Discharge from the slide occurs from Adi t 1 at about 4.2 x 10~ 2 m 3 / sec and 6.1 x 1 0 ~ 3 m 3 / s ec from Adi t 2. The remainder of the discharge occurs at the Revelstoke Reservoir and to channel deposits which are not continuous across the toe of the slide (Freeze, 1982). 4) Assumptions in the Analysis Figure 2.7a shows a plan view of the slide wi th boreholes completed since 1965 and some of the principal instrumentation locations. Bo th parts of Figure 2.7 show the location 2.12 Analysis of F ie ld Data of drainage adits and galleries excavated to reduce pore-pressures across the major-shear zone. Experience wi th the drainage system suggests that the transient response time of the slide is quite fast (i.e. steady-state behaviour is reached within weeks). Pa t ton (1984) describes the type of water level piezometers installed at the slide, and provides a review of the accuracies obtainable at some of the individual sites. These piezometers range from sealed tube-type standpipes to leaky inclinometer casing. When the details of all the instrumentation are examined, we find that the principal source of inaccuracy in the water levels is in the completion point of the piezometer. Thus, the elevation at which a particular hydraulic head measurement applies is uncertain. In some cases, deviations as large as ± 7 0 m (Patton, 1984). Table 2.6 lists the individual measurement point coordinates, hydraulic heads and uncertainties in measurement. Values of hydraulic head error are obtained from Pat ton (1982) or estimated from completion details at piezometer sites. For example, if the length of the completion zone is 10 m, then an uncertainty of ± 5 m is assigned to that measurement. These numbers are best viewed as relative weights between observations. In spite of these uncertainties, Downie Slide is considered one of the best instrumented rock slides in the world. Other assumptions in our analysis are: 1. Two-dimensional approximation of the flow system. In many geotechnical investiga-tions, the surface relief in the third dimension is considered small compared to the relief of the two dimensional section. The topography of the Downie Slide is typi-cal of this situation. Based on previous hydrogeologic studies (Freeze, 1977, 1982) groundwater flows are expected to be primari ly parallel in the third direction, thus a two-dimensional approximation seems val id . 2. Steady-state approximation of the flow field. The hydrogeology of the slide is subject to transient responses due to the seasonal nature of the influx to the system. A n examination of the piezometric records show that many of the larger fluctuations are confined to points close to the water table, and are small compared to the depth of 2.12 Analysis of Field Data the flow basin. These transients pass through the system within a few weeks. Ad i t discharges are relatively constant throughout the year, usually only varying within 10% of their base values. Therefore, it is reasonable to assume that the flow system can be conceptualized as being at steady state. 3. Projection of boreholes off-section. This topic is related to the first assumption, but deals wi th measurement uncertainties. Because of projections, measurement points wi l l not be at exactly the same elevation on different sections. This factor introduces an unknown component of noise into the measurements. These sources of noise are considered small compared to uncertainties in the completion points in the borehole. 4. Uncertain boundary conditions. Inverse techniques have most commonly been applied to areal aquifer models. In general, boundary conditions are more easily defined in these cases. On the particular section under consideration at Downie Slide, there are uncertainties in water table location on about one-third of the slide, and in this area we assume a free-surface boundary condition, the elevation of which is determined during the inverse. Good estimates exist on the water table location near the toe of the slide and on the recharge zone at the upper end of the Slide. Water balance studies provide ranges of recharge in the upper portion of the slide and continuous discharge records are available from both drainage adits. Assumed impermeable boundary conditions at the base of the model are believed to be far enough removed from the surface topography so as not to influence the flow system. A hydrostatic head distribution along C — C is assumed from regional flow considerations (Figure 2.6). Uncertainties in all boundary conditions (free-surface, input/output flux, etc.) wi l l be incorporated into the inverse. 5) Interpolation of Hydraulic Head Data Usually, an init ial step in applying the inverse procedure in aquifer problems is to estimate the values of hydraulic head at all nodal points in a finite element grid (for example, Cl i f ton and Neuman, 1982). The purpose of the interpolation is to augment the data base w i th both estimated values and a prior covariance structure for the hydraulic 2.12 Analysis of Field Data head data. This is a Bayesian approach to the inverse as noted in Appendix A . The inverse is then further conditioned compared to a non-interpolated data set by incorporating additional degrees of freedom. The difficulty in applying a linear estimator (as in kriging) to hydraulic head data is that a drift component exists in the data. In many examples in the literature the drift is successfully modeled by simple polynomial functions or removed by generalized increments (typical of universal kriging). To improve drift identification in more difficult problems, Neuman and Jacobson (1984) adopted a residual kriging scheme after representing the drift by a fourth- order polynomial . On Downie Slide, most of the data points are collected near the water table and are concentrated vertically. Futhermore, wi th the steep topography of this vertical section a strong drift component exists in the data, which cannot be adequately represented by a simple polynomial form. These two facts make the application of kriging difficult. In an attempt to deal with the limitations of the data set, we have adopted a more physically-based approach to kriging. A non-linear form for the drift is assumed and residual kriging is used. The form of the drift is given by the solution to Poisson's equation on the section, with two unknown coefficients that represent flux into the system along the upper boundary and flow out of a point sink. Adopt ing this form for the drift results in a better approximation to the average flow-field in the lower portions of the basin than is possible wi th a simple polynomial. The problem of identifying the two unknown parameters which determine the drift is non-linear and so the optimization procedure described earlier is used. Parameters are determined that yield a set of minimum residuals in hydraulic head. The algorithm is described in Appendix 2B. The data set suffers from large uncertainty in the measurements and inconsistency between measurement intervals in the vertical and horizontal directions. Because the flow is pr imari ly parallel to the mountain slope and the measurements are concentrated vertically, individual points may not contribute independent information to the drift determination. To test the ability of residual kriging to produce an unbiased estimate of the hydraulic 2.12 Analysis of Field Data head field, we perform a number of test runs on a synthetic data set. A hydraulic field is overlain on the cross-sectional geometry shown in Figure 2.8. Hydraulic conductivity values are generated for a lognormal distr ibution, wi th a mean log conductivity equal to -7.00, and a standard deviation of 0.315. A n exponential correlation function is adopted, wi th an integral scale of 200 m. It should be noted that the log-conductivity field is first generated on a regular grid, and then the values are transferred to the finite element mesh. In this case the vertical correlation structure of the original population is not preserved. However, this simplification does not affect the purpose for which the example is designed. The steady-state problem is solved by sub-dividing the problem area into 310 quadrilateral elements (620 triangles) wi th 352 nodes. The hydraulic heads at 33 points (in similar positions in the flow system as on Downie Slide ) are sampled and corrupted wi th Gaussian noise of N(0, 10 m). The residual kriging procedure is followed and the drift is subtracted from values at the sampled points. A variogram is computed for the residual data (Figure 2.9 ). The variogram quality is poor; obvious gaps in the data are present and low numbers of data pairs are indicated. The variogram appears to be uncorrelated with a sill of approximately 85 m 2 . When the data are kriged at all nodal points the kriging errors should be coherent with the predicted variance. In Table 2.7, E(h* — hact) is the average difference between the kriged and actual hydraulic heads. This number should be close to zero if the interpolation is unbiased. Sh* represents the noise level in the hydraulic head data set (equation 2.35). Sk is the average standard deviation in hydraulic head obtained from the kriging algorithm. Note that these two statistics should be equal if the kriging interpolator correctly determined the actual hydraulic head values. This group of statistics can be viewed in terms of two probabilty distribution functions: 1. A desired probabilty function with a mean of 0 and standard deviation of S A *. 2. A function produced from kriging with a non-zero mean and standard deviation of Sk-These two distributions should not be statistically different. If this criterion is not met, then any inverse based on the interpolated data wi l l produce biased parameters. Table 2.7 2.12 Analysis of Field Data (Trial 1) shows that there is a statistically significant difference between the two distribu-tions, indicat ing that the kriging procedure produces a biased estimate. The reasons for the poor performance of the interpolator can be attributed to the effects of noise in the data and to the large horizontal to vertical spacing of the data points. Wi th most of the data situated close to the water table, the estimate of the drift in the lower portions of the flow domain is inaccurate (de Marsi ly , 1986). To confirm these findings, a further simulation is carried out wi th the synthetic data. Here we sampled hydraulic head values at 33 data points distributed regularly in the X direction, wi th the depth of each position chosen randomly. The data is corrupted with noise, and the residual kriging operation was again followed. The resulting variogram (Figure 2.10) now indicates a linear covariance function with a sill of 850 m 2 . W i t h data points scattered throughout the domain, the variogram calculation now incorporates the larger variabil i ty from data pairs spaced further apart, leading to the larger value for the s i l l . When the data are kriged at all nodal points (Trial 2, Table 2.7) an unbiased estimate is obtained. Kr iged hydraulic heads and sample point locations are shown on Figure 2.11. It can be concluded from these simulations that a similar kriging exercise for the Downie Slide hydraulic data set would produce biased estimates. Hence, we cannot interpolate hydraulic heads to all nodal points and adopt a stochastic framework when solving the inverse problem. Therefore, the data base is limited to only 31 points (Table 2.6). The fact that no prior hydraulic conductivity data is available further limits the analysis. This restriction forces us to find the simplest model that fits the data. In the case 'simplest' refers to a homogeneous, isotropic hydraulic conductivity distribution. In the next section the inverse problem is solved for a hydraulic head data set collected at Downie Slide. To show the flexibilty of our inverse algorithm a stochastic inverse on the noise-free synthetic data set just described wi l l also be carried out. 6) Application of The Inverse Technique Based on the findings in the previous section, it is considered that it is not possible to 2.12 Analysis of Field Data generate an unbiased estimate of hydraulic head for Downie Slide using the existing data base. By examining both the experiments in the preceding section and Table 2.7, it is probable that local deviations in hydraulic head caused by spatial variations in hydraulic conduct ivi ty are about the same or smaller in magnitude than the noise level in the data. In applying the inverse method at this site, uncertain boundary conditions and uncertain geometries exist in the form of a free-surface over about one-third of the water table surface. Thus, the best that can be achieved in this case is to determine the simplest model that satisfies the data. For this problem the simplest model takes the form of a homogeneous, isotropic flow field, and the inverse problem takes the form of determining: the sink strength of the adit (modeled as a point), the free-surface location, a constant hydraulic head distribution along boundary C — C , and the fluid flux across the upper recharge area of the slide (Figure 2.12). This inverse could be classified as deterministic (weighted least squares). W i t h a free-surface present, it is not possible to use conventional inversion schemes to solve this type of problem. However, with the optimization technique proposed, the inverse posed on Figure 2.12 can be set-up and solved quite readily. The model and its unknown parameters are shown on Figure 2.12. The segment from A to C represents the inital estimate of the free-surface location, wi th the input flux parameter P\ applied from B to C. Along C - C values of hydraulic head are to be determined. Parameter P2 refers to the value of hydraulic head specified at the base of the flow system at C . A linear variation in hydraulic head from C ' to point C is then specified. This condition represents a simple form for this boundary, and is considered reasonable based on our examination of the regional flow system. Parameter P3 represents the sink strength of the adit system. Upper and lower bounds for all of the parameters are listed in Table 2.8. A n ini t ial estimate of the parameters is also listed, which is derived from forward model calibration. Note that the fluid fluxes have all been divided by an assumed hydraulic conductivity value of 1.0 x 1 0 - 7 m/sec. In fact, the inverse determines ratios of flux to hydraulic conductivity. In this simple model, the best-fitting hydraulic 2.12 Analysis of Field Data 67 conductivity value can only be determined if boundary fluxes are known with certainty. The functional & in this case is the solution to Poisson's equation with four model pa-rameters as described above. This equation is solved using the boundary element method ( B E M : Brebbia , 1978). This method is ideally suited for free-surface problems in homo-geneous media. A further advantage of the B E M is that it is only necessary to compute the potential at the specific measurement points and not throughout the entire domain as required by finite elements. This feature saves computational time when calculating functionals (2.7) or (2.11). The basic procedure as noted in section 2.10 is followed with the exception that forms of (2.7) and (2.11) without log-conductivity penalty terms are minimized. A n inner-iteration loop can be defined, in which the free surface is adjusted under each parameter set until the elevation of the free surface equals the hydraulic potential. The matrix Vh in each case is taken as a diagonal matrix with terms developed from Table 2.6. Note that this inverse problem is similar to the drift determination followed earlier in this chapter and detailed in Appendix 2B . However, here more parameters are determined, and the problem is posed for a different purpose. The results of the inversion runs are shown on Table 2.9. A total of 20 t r ia l points were used in the search procedure. Typical runs took 20 minutes on a F P S 1 6 4 / M A X array processor. Results presented in Table 2.9 are averages of the 20 parameter sets determined from the inverse. Standard deviations of the model parameters are also listed. In order to reduce the magnitudes of these values, more accurate hydraulic head data would have to available for the inverse. Recall that each of the trial sets produces data which satisfies the observations. Similar results for the two different norms indicates the strong possibility that these three final parameter estimates are good representors of the ' true' model. If all models had produced unacceptable misfits, for instance if the value of the objective function was not reduced below a level expected from the noise level in the data, then a satisfactory model would not have been found. Addi t iona l parameters would then have to 2.12 Analysis of Field Data be introduced, increasing the complexity of the model. A prudent modeler may ask: " If one possessed noise-free hydraulic head data at better locations wi thin the domain, how well could the hydraulic conductivities be determined using the approach suggested in this paper? " The question cannot be answered directly for the Downie Slide, but inferences in this regard can be made by working with synthetic data. In the preceding section an unbiased hydraulic head data set was generated when data points were located throughout the domain. Figure 2.11 shows the kriged hydraulic heads, and Figure 2.13 shows the estimation errors derived from kriging. Given the location of the water table, recharge along the upper surface, fluid discharges from the adit, and a known constant head boundary on the right hand side, the unknown log-conductivities in each of the 310 finite element blocks can be estimated. In two inverse simulations (Table 2.10) the L2 problem (2.7) and the L j problem (2.11) are solved with the penalty conductivity terms, using the hydraulic head estimation errors derived from kriging. Here, it is assumed that log-conductivity prior information is available in each block. The actual log-conductivity values are corrupted with Gaussian random noise with a standard deviation of ± . 2 0 . These values are then used as prior information in the inverse. The number of search points in this case is 400, the a priori range on the log-conductivity parameters is -9.00 to -5.00, and A is set equal to 1. Average runs took 25 minutes on the F P S array processor. Figures 2.14 and 2.15 show the h field determined from the Li inverse and the standard deviation in hest derived from cov hest. Note the similarity between Figures 2.11 and 2.14, and the reduction in the standard errors as a result of the inverse (Figures 2.15 and 2.13). The small values for the standard deviations in hydraulic head reflect the large data base used in this example. The log-conductivities determined by the Li inverse are shown on Figure 2.16. The L\ inverse produced identical results. The actual log-conductivity values are shown on Figure 2.17. For reasons of clarity, only two contour values; -7.2 and -6.8 are plotted on these two figures. The results show that the mean conductivity and the standard deviation of the original log-conductivity distr ibution are recovered (Table 2.10). 2.13 Discussion and Conclusions The standard deviation about the actual hydraulic conductivities (Table 2.10, Column 5) show that the computed log-conductivity field is determined to an accuracy approximately equal to that of the noise level in the a priori data. This result indicates that with known boundary conditions, and well-defined prior information, it is possible to determine unique parameterizations and aquifer statistics with the technique proposed. O f course, the preceding example uses a level of prior information not commonly en-countered in problems of this nature. In order to apply the method in a more realistic set-ting, another inverse is carried out; this time it is assumed that noise-free log-conductivity values are available only at the 31 measurement points as noted on Figure 2.11. In this simulation, the smallest deviatoric norm is minimized for a series of A values (equation 2.9). Vector Y" for these runs consists of 31 'observed' values and 279 terms set to a spec-ified value of -7.00. When minimizing (2.9), convergence is slower with larger values of A, but if the data misfit is satisfactory, parameter values replicate well between simulations with different random-number seeds. On the other hand, a value of A which is smaller leads to an algori thm which is faster, but which could yield non-unique results between different starting seeds. The best estimate of the Y field is shown on Figure 2.18. When Figure 2.18 and Figure 2.17 are compared, some of the major hydrogeologic features are replicated. Not all features are identified by the inverse; however, this is not surprising when considering that a deviatoric norm (equation 2.9) tends to smooth or damp varia-tions in hydraulic conductivity. This model (Figure 2.18) is a reasonable approximation to the ' true' model based on the data available. 2.13 Discussion and Conclusions In order to solve the inverse problem in groundwater hydrology two concepts need to be considered. Firs t , a formalism or philosophy about the problem needs to be developed. In this thesis, the philosophy is adopted that a variety of models should be constructed under different norms. Only the features of the model that are consistently replicated 2.13 Discussion and Conclusions are candidates for the 'true' model. Second, a robust inversion scheme capable of solv-ing the norms described above is needed. A groundwater inverse technique based on an opt imizat ion procedure described by Box (1965) and Silva and Hohmann (1983) is de-scribed. The algori thm requires no derivative computations and can be used to minimize an arbitrari ly-complicated non-linear functional. T h e algorithm produces a wide variety of acceptable models clustered about a global min imum. Each model generates data that matches observed values wi th in a noise level. Under the assumption that enough points have been sampled, the centroid of the sampled points is taken as the best estimate of the parameters. The covariances of the estimated parameters and estimated data are easily computed. A n y norm can be minimized with the algori thm. Various forms of both the L\ and L2 norms are chosen depending on assumptions of the nature of noise and a priori information available. It is customary in many inversion schemes to employ kriging techniques to interpolate hydraulic heads to all nodal points in a numerical g r id . Under noise-free conditions with adequate spacing between data points this interpolation is possible. It was observed how-ever, that on a cross-sectional problem with noisy data and a strong drift component, this interpolation was not possible. The difficulties are related to the inability of the residual kr iging scheme to correctly identify drift in the presence of noise with a poorly distributed data set. O n a practical problem, the technique has been applied to a sparse, inaccurate hy-draulic head data set at Downie Slide, Bri t ish Columbia . Results of a pre-inverse analysis indicate that the best that can be achieved at this site is to determine a simple model structure that explains the data. While more complicted structures may exist that explain the data, the simplest structure may be the most geologically unique. Four parameters are determined: the free-surface position of the water table and three boundary conditions for the domain. Further simulations using a theoretical data set with assumed properties s imilar to that of Downie Slide show that wi th noise free data, and an adequate spacing 2.13 Discussion and Conclusions between points it is possible to interpolate an unbiased estimate of hydraulic head data at all nodes in the equivalent discretized domain. When the inverse technique is applied, the domain's conductivity structure is correctly identified when enough prior log-conductivity information is available. The implications for Downie Slide are that in order to construct anything but a simple hydrogeologic model, accurate field measurements of hydraulic head are required, as well as well-defined estimates of hydraulic conductivity, a better spacing between measurements, and adequate knowledge of the boundary conditions. In Chapter 3, an additional data set (temperatures) wi l l be introduced in order to further condition the inverse problem. Further inverse simulations of Downie Slide wil l be carried out incorporating a temperature data set to assess possible benefits. Appendix 2A Appendix 2A T h e inverse approach is viewed as stochastic or deterministic depending on how the nature of the uncertainties in the data, the model, and the noise vector u are viewed in (2.2). If the noise vector v or the model is viewed as random, then the inverse is stochastic. For other cases, the model is deterministic. Most of control theory (e.g. Schweppe, 1973) categorizes these inverses according to the level and nature of a priori information available. These categories are: (1) Bayesian: The model p, and the vector v are viewed as random vectors. Mathemat i -cal ly we can express this as: E{p) = Pz E\{P - Pz){P - Px)T] = COV P" = OyVy E(u) = 0 (Al) E(vuT) — cov h' — o\Vh E(u,p) = 0 If (2.7) is minimized then the model takes the form of unknown Y values, and cov p~ = <7yVy•, which can be estimated by kriging. One may also adopt an assumed form for V y and minimize (2.7) for a range of A values (Gavalas, et. a l , 1976). (2) Fisher: The vector p has no a priori statistical structure, and v is treated as a random vector. In this case: E(u) = 0 (A2) E(vvT) — cov h' (3) Weighted Least Squares: Both the v and p are completely unknown a pr ior i . Generally a simple form of (2.7) is minimized (i.e.): J = (h — h*)TW(h — h*) (A3) Appendix 2A were W is a positive-definite weighting matrix; the terms of which are usually chosen by judgement or empirical rules. (4) Unknown but bounded: The model p and the noise vector v are unknown but are constrained to lie wi thin certain bounds. A n example of this type of estimator would be the minimizat ion of functional (2.10), with upper and lower limits specified in advance. Appendix 2B Appendix 2B The development of the residual kriging scheme is similar to that of Neuman and Jacobson (1984). Let Z be a vector of noise-free Z(x) values measured at / discrete locations, xl,i = 1,2,3 and let n be the corresponding drift vector. Z and ^ are related by •• Z = fi+'R (51) where R is a vector of residuals at xt. We can characterize R by a covariance matrix V , defined as E\RRT\ = V (B2) The scheme is iterative in nature, and let p represent the current iteration level. The drift at the pth step of the procedure is represented as: 2 lxp{x) = Y,H{x,aj) ( £ 3 ) j = i where H(x,a,j) is a nonlinear function given by the solution to the homogeneous, isotropic flow equation V2H = aiS{x- x') (54) and the two unknown parameters are a\ the sink strength, and a 2 the flux into the system along the top boundary. If V is known, we can obtain an unbiased estimate of the aj by minimizing n{a) — {Z - p,)TV~l(Z - p.) (B5) where d are the values of the parameters computed from the solution to (B5) and p, = H(x,a). The solution to (B5) can be found by using the constrained simplex procedure described previously. The kriging procedure is as follows: A t p = 1, V is chosen as a diagonal matrix, and the jl drift function is computed yielding an estimate R of the residuals. If the variogram Appendix 2B of the residuals ^(R) has a distinct sill and no anisotropy, then we proceed to kriging the residuals (Neuman and Jacobson, 1984). If no sill is present we estimate a new V matrix by estimating a sill and where 7^ is the estimated si l l , s tj = \xl - Xj\ , and p^(0) is the covariance at zero lag. The p variable is incremented by one and (A4) is minimized again. Extensions to kriging noisy data can also be implemented (Journel and Huijbregts, 1978). C h a p t e r 2 - N o t a t i o n 76 Notation A global stiffness matrix B new position of point in simplex figure cov hesi estimated hydraulic head covariance matr ix cov h~ a priori hydraluic head covariance matrix cov mest estimated model covariance from inverse d vector of computed data d~ vector interpolated or observed data values 5 functional transformation / vector of boundary conditions C7,y sensitivity matr ix G~9 generalized inverse of G g(x,y,z) function specifying Dirichlet conditions on T 2 H non-linear function of drift H*k Storage matrix: k ' th estimate of i ' th value of hydraulic head h hydraulic heads computed by finite element. h~ observed or interpolated (from kriging) hydraulic heads hest estimated hydraulic heads from inverse, based on mest h centroid of hydraulic heads from inverse average value of the i ' th hydraulic head I identity matrix / number of data points J objective function K hydraulic conductivity tensor L lower triangle matr ix in Cholesky scheme L number of search points L\, L2 vector norms referring to functional types Chapter 2 - Notation 77 M number of observed and estimated data points m N dimensional vector of parameters ./V number of parameters Po m in imum length solution p- 1 arbitrary vector of null-space components P^k k ' th estimate of i ' th parameter p~ average value of the i ' th parameter Q previous position of point in complex q sources or sinks in the groundwater flow equation qj specified flux on the boundary q' specified dimensionless fluid-flux on the boundary R vector of residual heads S centroid of the remaining points in the simplex S~h root mean square deviations from inverse to actual data S/i- root mean square deviations from observed to actual data Sy root mean square deviations from inverse to actual value Sy root mean square deviations from observed to actual value V scaled covariance matrix of residual heads Vh scaled covariance matr ix of estimated hydraulic heads Vy scaled covariance matrix of estimated log conductivity x vector of grid coordinates in a numerical scheme Y, Y, Y* log-hydraulic conductivity, centroid of all Y values from the inverse, and observed or interpolated (from kriging) log-conductivity Y / , Yu prior constraints on log conductivity in the form of lower and upper bounds Z vector of observed data in the residual kriging scheme a expansion factor in constrained simplex scheme Ti, T 2 boundary domain segments C h a p t e r 2 - N o t a t i o n 7 integral scale of a random variable l{sij) semi-variogram at lag s A scaling parameter in multi-objective criteria (j, vector of drift or trend of h at a set of disecrete points v vector of residuals between computed and observed hydraulic head \I> theoretical hydraulic head Pfi(0) covariance at zero lag ah scaling factor used in defining magnitudes of covariances in hydraulic head ay scaling factor used in defining magnitudes of covariances in log conductivity ox2 standard deviation of the chi-squared variable X2 chi-squared statistic X 1 chi-one statistic Chapter 2 - References References Box, M . J . , 1965. A new method of constrained optimizat ion and a comparison wi th other methods, Computer Journal, 8, 42-52. Brebbia , C . A . , 1978. The Boundary Element Method For Engineers. Pentech Press. Bir t les , A . B . and E . H . Morel , 1979. Calculat ion of aquifer parameters from sparse data, Water Resources Research, 15(4), 832-844. Carrera , J . and S.P. Neuman, 1986a. Est imat ion of aquifer parameters under transient and steady-state conditions, 1, M a x i m u m likelihood method incorporating prior infor-mation, Water Resources Research, 22(2), 199-210. Carrera, J . and S.P. Neuman, 1986b. Est imat ion of aquifer parameters under transient and steady-state conditions, 2, Uniqueness, stability and solution algorithms, Water Resources Research, 22(2), 211-227. Carrera, J . and S.P. Neuman, 1986c. Est imation of aquifer parameters under transient and steady-state conditions, 3, Appl icat ion to synthetic and field data, Water Resources Research, 22(2), 228-242. Cl i f ton, P . M . and S.P. Neuman, 1982. Effects of kriging and inverse modeling on condi-tional simulation of the A v r a Valley Aquifer in Southern Ar izona , Water Resources Research, 18(4), 1215-1234. Cooley, R . L . , 1977. A method of estimating parameters and assessing reliability of models of steady state groundwater flow, 1. Theory and numerical properties, Water Resources Research, 13(2), 318-324. 1 Cooley, R . L . , 1979. A method of estimating parameters and assessing reliability of mod-els of steady state groundwater flow, 2. Appl ica t ion of statistical analysis, Water Resources Research, 15(3), 603-617. Cooley, R . L . , 1982. Incorporation of prior information on parameters into nonlinear regres-sion groundwater flow models, 1. Theory, Water Resources Research, 18(4), 965-976. Cooley, R . L . , 1983. Incorporation of prior information on parameters into nonlinear re-gression groundwater flow models, 2. Applicat ions, Water Resources Research, 19(3), 662-676. Cooley, R . L . and R . L . Naff, 1985. Regression modeling of groundwater flow. U S G S Open File Report 85-180. Dagan, G . , 1985. Stochastic modeling of groundwater flow by unconditional and condi-tional probabilities: The inverse problem, Water Resources Research, 21(1), 65-72. Dal ton , R . D . , 1986. Application of a new finite element program to the analysis of ground-water flow in Downie Slide, British Columbia. B . A . S c . Thesis, Department of Geo-logical Sciences, University of Bri t ish Columbia . C h a p t e r 2 - Refe rences Draper, N . R . and H . Smith , 1981. Applied Regression Analysis. Wiley & Sons, New York. Freeze, R . A . , 1977. Computer modeling of groundwater conditions on Downie Slide. Re-port to Hydroelectric Design Divis ion , Br i t i sh Columbia Hydro, Vancouver, Canada. Freeze, R . A . , 1982. Groundwater conditions on Downie Slide, 1, Review of the available hydrogeological data on Downie Slide. Report to Hydroelectric Design Div i s ion , Bri t ish Co lumbia Hydro , Vancouver, Canada. Freeze, R . A . and J . A . Cherry, 1979. Groundwater. Prentice-Hall . Fullagar, P . K . and D . W . Oldenburg, 1984. Inversion of horizontal loop electromagnetic frequency soundings, Geophysics, 4 9 ( 2 ) , 150-164. Gavalas, G . R . , P . C . Shah and J . H . Seinfeld, 1976. Reservoir history matching by Bayesian estimation, Soc. Pet. Eng. J., 16, 337-350. Hoeksema, R . J . and P . K . Ki tan id is , 1985. Comparison of Gaussian conditional mean and kr iging estimation in the geostatistical solution to the inverse problem, Water Resources Research, 21(6) , 825-836. Imrie, A . S . and D . R . Bourne, 1981. Engineering geology of the M i c a and Revelstoke Dams, Field guides to geology and mineral deposits. In Proc. Joint Meeting GAC, MAC, and CGU, (393-401), Calgary. Journel , A . G . and C . J . Huijbregts, 1978. Mining Geostatistics. Academic Press, New York. Ki tan id i s , P . K . and E . G . Vomvoris, 1983. A geostatistical approach to the inverse prob-lem in groundwater modeling (steady state and one dimensional simulations), Water Resources Research, 19(3) , 677-690. de Mars i ly , G . , G . Lavedan, M . Boucher, and G . Fasanino, 1984. Interpretation of in-terference tests in a well field using geostatistical techniques to fit the permeability dis tr ibut ion in a reservoir model. In Geostatistics For Natural Resource Characteriza-tion, Part 2, (G . Verly et al. , eds.), Reidel . de Mars i ly , G . , 1986. Quantitative Hydrogeology. Academic Press. Meneley, W . A . , 1977. Interpretation of the recovery test on drillholes U-5 and U-6, Downie Slide. Repor t to Hydroelectric Design Divis ion, Br i t i sh Columbia Hydro, Vancouver, Canada. Menke, W . , 1984. Geophysical Data Analysis: Discrete Inverse Theory. Academic Press, New York . Neuman, S.P., 1980. A statistical approach to the inverse problem in aquifer hydrology, 3. Improved solution technique and added perspective, Water Resources Research, 16(2) , 331-346. Neuman, S.P., 1982. Statistical characterization of aquifer heterogeneities: A n overview. In Recent Trends in Hydrogeology, ( T . N . Narasimhan, ed.), G .S .A . Special Paper 189. Chapter 2 - References Neuman, S.P. and S. Yakowitz, 1979. A statistical approach to the inverse problem in aquifer hydrology, 1. Theory, Water Resources Research, 15(4), 845-860. Neuman, S.P., G . E . Fogg and E . A Jackobson, 1980. A statistical approach to the inverse problem in aquifer hydology, 2. Case study, Water Resources Research, 16(1), 33-58. Neuman, S.P. and E . A . Jacobson, 1984. Analysis of nonintrinsic spatial variability by residual kriging with application to regional groundwater levels, Math. Geol, 16(5), 499-521. Oldenburg, D . W . , 1984. A n introduction to linear inverse theory, IEEE Trans. Geosc. and Remote Sens, GE-22(6), 665-674. Parker , R . L . , 1977. Understanding inverse theory. A n n . Rev. Ear th Planet. Sci..35-64 Parker, R . L . and M . K . M c N u t t , 1980. Statistics for the one-norm misfit measure, J. Geophys. Res., 85(B8), 4429-4430. Pa t ton , F . D . , 1983. The role of instrumentation in the analysis of the stability of rock slopes. In Proc. Int. Sympos. on Field Measurements, in Geomech., Vol 1, (719-748), Zur ich . Pa t ton , F . D . , 1984. Climate, groundwater pressures and the stability analysis of landslides. In Proc. Int. Sympos. on Landslides, Vol. 8, (1-17), Zurich. P i teau , D . R . , F . H . Myl rea and I.G. Blown, 1978. Downie Slide. Bri t ish Columbia, Canada. In Rock Slides and Avalanches, Vol. 1, ( B . Voight. ed.), 365-395. Pr ice , W . L . , 1965. A controlled random search procedure for global optimization, Com-puter Jour., 8, 42-52. Reklai t i s , G . V . , A Ravindran and K . M . Ragsdell , 1983. Engineering Optimization, Meth-ods and Applications. Wiley, New York. Schlax, M . G . , 1984. A practical solution for linear inference computations. M . S c . Thesis, Department of Geophysics and Astronomy, University of Brit i tsh Columbia. Schwefel, H . -P . , 1981. Numerical Optimization of Computer Models. Wiley. Schweppe, F . C . , 1973. Uncertain Dynamic Systems. Prentice-Hall. S i lva , J . B . C . and G . W . Hohmann, 1983. Nonlinear magnetic inversion using a random search method, Geophysics, 48(12), 1645-1658. Tarantola, A . , and B . Valette, 1982. Generalized nonlinear inverse problems using the least square criterion, Ann. Rev. Geophys. Space Phys., 20(2), 219-232. Townley, L . R . and J . L . Wilson, 1983. Condi t ional second moment analysis of groundwater flow: The cumulative effects of transmissivity and head measurements. In Proc. Int. Conf. Groundwater and Man, (Sydney), Austral ia . Wiggins , R . A . , 1972. The general linear inverse problem: Implication of surface wave and free oscillations for earth structure, Rev. Geophys. and Space Phys., 10(1), 251-285. Chapter 2 - References Yeh, W . W - G . , 1986. Review of parameter indentification procedures in groundwater hy-drology: The inverse problem, Water Resources Research, 22(2), 95-108. Yeh , W . W - G . , Y . S . Y o o n , and K . S . Lee, 1983. Aquifer parameter identification with kriging and opt imum parameterization, Water Resources Research, 19(1) , 225-233.3. Chapter 2 - Tables Table 2.1 Example construction problem - Four zone model with noise free data Parameter L\ L2 Actual Pi 10.0 10.0 10.0 P2 1.0 1.0 1.0 P 3 1000. 1000. 1000. Data at measurement points: hi = 137.62 h2 = 137.55 h3 = 137.48 h4 = 136.81 h5 = 136.13 h6 = 136.07 h7 = 136.00 / i 8 = 68.51 hg = 1.00 Chapter 2 - Tables Table 2.2 Example construction problem - Four zone model wi th blunder in eighth data point Parameter L i L 2 Actual P i 10.0 384.05 10.0 P2 1.0 -181.6 1.0 P3 1000. 1156.3 1000. D a t a at measurement points: hx = 137.62 h2 = 137.55 h3 = 137.48 h4 = 136.81 h5 = 136.13 hG = 136.07 h7 = 136.00 hg = 130.00 hg = 1.00 Chapter 2 - Tables Table 2.3 Example construction problem - Two zone model wi th noise free data (Incorrect parameterization) Parameter L\ L2 Ac tua l P i 505.0 545.75 10.0 P2 - - 1.0 P3 - - 1000. Data at measurement points: hx = 137.62 h2 = 137.55 hz = 137.48 h4 = 136.81 h5 = 136.13 h6 = 136.07 h7 = 136.00 hs = 68.51 h9 = 1.00 Chapter 2 - Tables Table 2.4 One-dimensional, steady-state hydraulic-head inverse. Synthetic example wi th 20 nodes, 19 unknowns. L 2 norm minimizat ion, h" noise level in hydraulic head data in % of actual value, Y" noise level in log-conductivity data in % of actual value. ' N O ' is defined as having no prior information incorporated into the inverse. root mean square inverse-misfit, Sh- noise level in hydraulic head data, S y root mean square inverse-misfit, Sy noise level in prior log conductivity data. roup h" Y" Sh Sy A 10 N O .014 .074 .424 _ 10 10 .007 .078 .152 .202 10 20 .023 .074 .293 .400 B 10 N O .044 .063 .221 -10 10 .030 .063 .207 .210 10 20 .038 .063 .249 .412 C 10 N O .018 .075 .346 -10 10 .015 .075 .114 .197 10 20 .018 .075 .251 .400 Group A : aY = 1.0, 7 = 0.2, Group B : aY = 0.1, 7 = 0.2, Group C : aY = 1.0, 7 = 0.8. Chapter 2 - Tables Table 2.5 Comparison between Carrera and Neuman (1986c, Table 2.1) and the constrained simplex technique (CT). Theoretical example with no prior information. Hydraulic heads corrupted with N (0, 1 m) noise. Transmissivity values in m 2/day. Zone Actual C & N LogError C T 1 LogError C T 2 LogError 1 150 125.74 0.077 156.49 -0.019 158.49 -0.024 2 150 201.70 -0.129 151.62 -0.005 158.49 -0.024 3 50 36.03 0.142 55.45 -0.044 44.24 0.054 4 150 156.45 -0.018 152.52 -0.007 133.05 0.052 5 50 42.26 0.073 52.83 -0.023 53.58 -0.029 6 15 17.26 -0.061 13.82 0.036 14.84 0.005 7 50 53.00 -0.025 47.15 0.027 47.87 0.020 8 15 14.74 -0.008 13.98 0.031 15.56 -0.016 9 5 4.87 0.011 5.29 -0.023 4.88 0.012 Sy .077 .027 .031 Note: Sy for Carrera and Neuman's other two realizations were 0.141 and 0.050. C T 1 , 2 refer to two separate realizations of noise added to the hydraulic head data. Chapter 2 - Tables Table 2.6 Downie Slide hydraulic head data set Standard errors o^- estimated from piezometer completion details X ( m ) Y ( m ) Head(m) ah.(m) 265 518 582 6 404 587 587 3 480 434 607 7 480 473 584 4 480 516 587 6 480 547 591 5 480 588 595 52 534 534 591 2 640 508 595 1 640 537 587 1 640 566 587 1 716 518 601 1 716 549 588 1 991 618 700 15 1381 907 951 76 1776 860 1070 12 1776 899 1142 70 1776 912 1145 45 1776 927 1114 46 1982 921 1143 15 1982 1025 1166 5 1982 1041 1177 6 1982 1072 1183 5 1982 1091 1183 8 1982 1118 1140 11 1982 1141 1152 23 2195 1111 1180 76 2713 1148 1379 42 2713 1160 1379 21 2713 1299 1308 8 2713 1315 1317 34 Chapter 2 - Tables 89 Table 2.7 Compar ison of kriging results for two different samples. Tr ia l 1 refers data gathered in a predominently vertical direction. Tr ia l 2 refers to data sampled regularly in the horizon-tal and random vertically. Sk is the average standard deviation (estimation errors from kriging) . S/,- noise level in hydraulic head data, Tr ia l E(h' - hact) Sh. Sk 1 2 4.43 -0.05 23.3 11.3 9.22 8.47 Chapter 2 - Tables 90 Table 2.8 Input parameter ranges for Downie Slide inversion Parameter Lower Bound Initial Estimate Upper Bound Pi Input fluid flux .050 .097 .120 P 2 Sink strength (m) -300 -240 -200 P3 Specified head(m) 1320 1331 1360 Chapter 2 - Tables Table 2.9 Comparison of results for the Downie Slide inverse Norm . Pi opv P 2 ^p 2 3^ op3 Li .075 .021 -218 18.3 1334 10 L2 065 .023 -220 18.1 1333 11 Note: Pi input fluid flux (m/m), P 2 sink strength (m), P3 specified head (m). Chapter 2 - Tables 92 Table 2.10 Synthetic Hydraulic Head Data Set Statistics of log-conductivity realization: //y = -7.00,oy = 0.315, S -^ = 7.63 Norm Sh E{Y) a y SY Li 5.06 -7.00 0.369 0.193 L 2 4.59 -6.99 0.346 0.181 Chapter 2 - Figures 10 Y 0 H r " — " — i 1 1 " l 0 2 4 6 8 10 X Figure 2.1 Graph of the curve y = 1.0 + 0.5i 2 , and four data points located along the curve. Both the L\ (dashed line) and L2 (solid line) solutions to the inverse are shown, y — 1.0+ 1.47x and y = 0.3 -f 1.96x, respectively. Chapter 2 - Figures h 9 - h L a. True model b. Incorrect model Figure 2.2 One dimensional groundwater flow system. Nine measurements of hydraulic head. Speci-fied fluid flux at z = 0 and h(L) = hi. (a) Actual model, consisting of four zones K\, K2, Kz, K4. Three parameter linear inverse, (b) Incorrect interpretational model consisting of two zones. Chapter 2 - Figures b. " m ••t model parameter m c . j i m M t mwt m"» m ; > t m M t m M t model parameter m model parameter m model parameter m nv m. Figure 2.3 (Adapted from Menke, 1984) Objective function J as a function of model parameters m. (a) Single minimum corresponds to an inverse problem with a unique solution, (b) Two solutions, (c) Many well seperated solutions, (d) Finite range of solutions, (e) Contour plot of J. There is an infinite range of solutions that lie anywhere in valley. Chapter 2 - Figures 1.0 r o.i o.e Q 0.4 I 0.2 0.0 -4.0 y • 0.2 fly " -2.0 Oy - 0.1 0>) -3.0 -20 -1.0 Log Conductivity o.o 1.0 0.0 8 >- 0.0 1 8 2 0.4 I N 0.2 0.0 y - o.80 fly - -2-0 (Jy " 1.0 J I I • ' -4.0 -3.0 -2.0 -1.0 Log Conductivity 0.0 (C) 0.8 E o.e e St 8 § 0.4 0.2 1± y - 0.20 fly - "2.0 Oy " 1.0 0.0 I I I I 1 L -4.0 -3.0 -2.0 -1.0 0.0 Log Conductivity Figure 2.4 Log-conductivity realization for one-dimensional examples: (a) Group A , (b) Group B , (c) Group C . Chapter 2 - Figures No Flow / y / / / / / / / / y y y / / / / / / / / / / / / / AREAL RECHARGE (I) 0.274 10^ 3 m/dov AREAL RECHARGE (2) 0.137 IQ-3 m/dav c 6 km PRESCRIBED HEAP (100m) \ h r r r r a. Flow Region 50 15 5 (7) (8) (9) 150 50 15 (4) (5) (6) 150 150 50 (1) (2) (3) b. Transmissivity Zones x 250(1) 250(1) i 250 (I) .-1000(3) . 250 (I).-1000(2) 1 250(1) \ 250 (I) (m3/day) c. Prescribed Flow Terms 16 17 18 1 5 .12 .13 " 4 •10 II • 9 .7 .8 • 4 .5 • 2 .3 .6 d. Observation Points Figure 2.5 Example two-dimensional inverse problem. Adapted from Carrera and Neuman (1986c). Chapter 2 - Figures 98 Figure 2.6 Regional groundwater flow system. Mode l bounded by Revelstoke Reservoir on the left (East) and the Monashee Mountains on the right (West). Chapter 2 - Figures metres Figure 2.7 Site map of Downie Slide, British Columbia, Canada, (a) Plan (b) Section. Chapter 2 - Figures 100 Figure 2.8 Finite element mesh: 352 nodes, 610 triangular elements. Chapter 2 - Figures 101 140 n 120 H 100 H 20 H 0 I i 1 1 r 0 200 400 600 800 Lag (m) Figure 2.9 Variogram for vertically oriented measurement points with noise. Chapter 2 - Figures 102 Lag h (m) Figure 2.10 Variogram for noise-free scattered data. Chapter 2 - Figures 103 1250 -i 0 500 1000 1500 2000 2500 3000 DISTANCE, metres Figure 2.11 Kriged theoretical hydraulic head data. Scattered data point locations. Chapter 2 - Figures 104 Figure 2.12 Downie Slide inverse problem: boundary conditions and parameters. P i : fluid recharge rate between points B and C. P 2 : constant hydraulic head at point C. P 3 : sink strength from adit system (modeled as point sink). Chapter 2 - Figures 105 Figure 2.13 Est imat ion errors in hydraulic head from kriging. Chapter 2 - Figures 106 1280 -i DISTANCE, metres Figure 2.14 Z/2 inverse hydraulic heads h. Chapter 2 - Figures 107 1280 -i DISTANCE, metres Figure 2.15 Li inverse standard deviations in h determined from cov h . Chapter 2 - Figures 108 1280 -i DISTANCE, metres Figure 2.16 Log-conductivi ty field from L2 inverse wi th prior information in each block. Non-shaded areas: log-conductivity between -7.2 and -6.8. Shaded areas: greater than -6.8. Lined areas: less than -7.2. Chapter 2 - Figures 109 1250 - i 1000 H DISTANCE, metres Figure 2.17 Actual log-conductivity field. Non-shaded areas: log-conductivity between -7.2 and -6.8. Shaded areas: greater than -6.8. Lined areas: less than -7.2. Chapter 2 - Figures 110 1280 -i DISTANCE, metres Figure 2.18 Smallest deviatoric log-conductivity field. Non-shaded areas: log-conductivity between -7.2 and -6.8. Shaded areas: greater than -6.8. Lined areas: less than -7.2. C H A P T E R 3 Hydraulic and Thermal Inversion S ummary A simultaneous inversion scheme involving both hydraulic, head and subsurface tempera-ture data is introduced. The joint inversion method uses data from a number of different data sets to improve the resolution of jointly held parameters. One such data set is sub-surface temperature, which is sensitive to variations in hydraulic conductivity in more permeable systems. The thermal and hydraulic head fields are linked by the fluid specific discharge. The basic idea behind the joint inversion scheme is simple; temperature mea-surements are taken in piezometers or boreholes along with hydraulic head measurements. M a n y more temperatures than hydraulic head measurements can be taken in most situ-ations. In addition, temperature measuring equipment is simple to construct and use in the field. The simultaneous inverse scheme is tested on a variety of synthetic examples. It can be concluded from studies of simple hydraulic conductivity structures that the joint inverse scheme is sensitive to the overall Bu lk Peclet number of the system, which is a measure of the level of advective disturbance. In purely conductive fields, inversion of the temperature field to construct hydraulic conductivity models does not yield good results, the temperature field being far too insensitive to hydraulic conductivity varia-tions. Similarly, at high fluid velocities variations in hydraulic conductivity does not effect temperature patterns. In other examples the ability of the simultaneous inverse scheme to reconstruct a hypothetical layered medium is examined. W i t h a data base limited to 111 Summary 112 hydraulic heads only, an inverse may indicate a smearing of the hydraulic conductivity in zones below measurement points and fail to determine boundary conditions such as the recharge rate to the flow system. In the synthetic examples chosen, the inclusion of thermal data allowed for the determination of the recharge rate to the flow system, and the location and hydraulic conductivity of an aquifer within the flow system. Addi t ional runs on a stochastic problem initiated in Chapter 2 were carried out. W i t h a temperature data base, thermal properties can be properly resolved. However, for the stochastic prob-lem (with known boundary conditions) the addition of thermal data did not condition the inverse to a greater degree than the addit ion of model prior information. The benefits of including thermal data and applying a joint inversion can be substantial when considering the more realistic problem of uncertain boundary conditions. The simultaneous inverse is also applied to the Downie Slide data set examined in Chapter 2. Results from that chapter show a simple model of a homogenous, isotropic medium. Four parameters were solved for: the free-surface position of the water table, and three boundary conditions for the domain. Unfortunately, with a homogeneous hydraulic conductivity, all that can be determined from the inverse are ratios of flux to hydraulic conductivity. B y including thermal data, the value of hydraulic conductivity can be determined at this site. Some of the model parameters (basal heat flux, thermal conductivity, specified head boundaries) are not resolved well by the joint scheme. None the less the constructed models do offer valuable insight into the hydrogeology of the field site. The constructed models persistently show a hydraulic conductivity structure of about 1 x 1 0 - 7 m/sec, which is consistent wi th previous estimates of hydraulic conductivity at the site. 3.1 Introduction 113 3.1 Introduction Difficulties associated with direct measurement of the hydrologic parameters needed for physically-based mathematical models are well known. Equally well known are the difficulties in the calibration procedure when trying to adjust parameters within a priori l imits unt i l model output at selected points matches observed values. Quite often questions are raised as to the uniqueness and optimality of these models. A major focus of research has been directed towards inversion techniques and parameter estimation as a way of both automatic calibration and as a statistical procedure to quantify the reliability of parameter estimates. In the inverse approach measurements of hydraulic head, hydraulic conductivi ty or transmissivity, seepage flux, and the like are inputs to the inverse algorithm, and fitted hydraulic conductivity (or other parameters) become the output, along wi th the parameter covariance structure. Throughout this paper the term 'model ' is defined as a configuration of material properties. Mathematically, the model takes the form of a vector p which consists of all unknown parameters (for instance, hydraulic conductivities, boundary conditions). Three separate aspects of an inverse problem need to be considered. First , specific goals about what is to be achieved from the inverse need stating. For instance, do we want the inverse to generate: (a) a set of models that reproduce the observed data, (b) unique averages of model parameters, or (c) upper and lower bounds to the model? Second, given the goal of an inverse problem, a suitable solution technique must be developed. The requirements of the technique should include: flexibility, robustness, and a minimum number of statistical assumptions either implici t ly or explicity implied in the method. T h i r d , consideration must be given to the nature of the observed data set. If the data are inaccurate and sparse, then model parameters may be poorly resolved. In the second chapter, a review of current methodologies of stating and solving the inverse was carried out, followed by the adoption of a formalism that advocates the con-struction of a variety of models under different norms. The inverse is carried out using 3.1 Introduction 114 a constrained simplex procedure as described by Box (1965). The algorithm requires no derivative computations and can be used to minimize an arbitrarily-complicated non-linear functional, subject to either linear or non-linear inequality constraints. The algorithm produces a variety of acceptable models clustered about a global min imum, each of which generates data that matches observed values. Under the assumption that enough points have been sampled, the centroid of these points is taken as the best estimate of the pa-rameters. In Chapter 2 this inverse approach is applied to the solution of steady state hydrogeologic inverses; including a data set from Downie Slide, Br i t i sh Columbia . In formulating a groundwater flow model, one can face uncertainties in assigning a spatial representation for the hydraulic conductivity field, K , identifying boundary condi-tions, and specifying the magnitude of fluid sources and sinks. One way of overcoming data l imitations in an inverse problem is to introduce a data set for a second potential field and apply simultaneous or joint inversion. The joint inversion method uses data from a number of different measurements to improve the resolution of parameters which are in common to one or more functional relationship. One such data set is subsurface temperature, which is sensitive to variations in hydraulic conductivity in more permeable systems (Smith and Chapman , 1983). The thermal and hydraulic head fields are linked by the fluid specific discharge. The idea of examining thermal fields to gain a better understanding of the hydrologic regime is best illustrated with a simple example. Figure 3.1 shows the thermal regime in a hypothetical basin for three different values of homogeneous, isotropic permeability: 1.0 x 1 0 ~ 1 8 m 2 , 2.0 x 1 0 - 1 6 m 2 , and 5.0 x 1 0 ~ 1 6 m 2 . The basin is 40 k m wide and 5 k m deep, wi th a linear water table having a total relief of 500 m . Heat flow within the earth is predominantly vertical in conductive environments. However, because of a coupl ing between the groundwater and heat flow equations, the thermal field wi l l show increased advective disturbance with increased fluid velocities. The upper plot in Figure 3.1 illustrates one end member in which the specific discharge is too low to have an effect 3.1 Introduction 115 on the conductive heat flow regime. The temperature within the basin in this simulation is governed entirely by the geometry of the basin, the thermal conductivity of the saturated porous medium, and the basal heat flux. As permeability is increased, fluid velocities become sufficient to redistribute heat in the system. Isotherms in the recharge area are depressed because of the downward flow of cooler water from the water table. In the middle of the basin, isotherms are t i l ted with respect to their conductive configuration, but for a broad region they remain subparallel to the surface. In this region, equipotentials are near vertical and have little sensitivity to variations in permeability. Isotherms in the discharge area are elevated because of the upward flow of warmer water at depth in the basin. A number of studies have demonstrated the sensitivity of the thermal field to vari-ations in the magnitude of hydraulic conductivity and anisotropic ratios, length/depth ratios of the flow system, thermal conductivity changes, varying basal heat flux, and the existence and location of aquifers (for example, Smith and Chapman, 1983). Three di-mensional effects of fluid flow on the thermal regime have also been examined (Woodbury and Smi th , 1985). Smith and Chapman (1983) describe an advective threshold, which is a relatively abrupt transition from a conductive to advectively-dominated system as fluid velocity increases with permeability. This transition occurs over approximately one order of magnitude in permeabilty for typical large-scale sedimentary basins. A t high fluid ve-locities the thermal regime of a basin can be dominated by advective effects, with a near isothermal temperature field. The basic idea behind the joint inversion scheme is to exploit the sensitivity of the ther-mal field to hydrogeologic parameters. Temperature measurements are made in piezome-ters or boreholes along wi th hydraulic head measurements. Normally only a l imited number of hydraulic head measurements can be taken in an individual piezometer. Temperatures, in contrast, can be profiled continuously along the length of a borehole or piezometer standpipe. This feature greatly increases the amount of temperature data compared to hydraulic heads that can be gathered at most sites. In addition, temperature measuring 3.2 Literature Review 116 equipment is simple to construct and use in the field. A n important l imitation st i l l remains; namely, limits on the number of boreholes in a basin. Temperature measurements are best made in fluid-filled boreholes that are in thermal equi l ibr ium with the rock mass. However, problems wi th temperature data sets can arise when fluid flow (either in or around a borehole) disturbs the thermal field at a measurement point. For this reason, boreholes should be cased or grouted prior to thermal profiling. These problems wi l l be discussed in more detail later in the Chapter. To estimate parameters for a groundwater flow model, a joint inversion scheme which minimizes a multi-objective criterion is carried out. Addi t ional unknown parameters such as thermal conductivity and thermal boundary conditions are also identified in the scheme. It is shown that for certain ranges of fluid flow, limitations in hydraulic head data sets, such as incompleteness or inaccuracy, can be overcome with the joint inverse scheme. In the following sections the basic concepts of heat and fluid transfer in porous media are reviewed, followed by the introduction of a series of objective functions and norms used for model construction. This introduction is followed by inversion of theoretical data and a re-investigation of the hydrogeology of Downie Slide, augmented with thermal data and a simultaneous inverse. 3.2 Literature Review A review of studies which consider various aspects of coupled fluid flow and heat trans-fer in regions wi th 'normal ' geothermal gradients is given by Smith and Chapman (1983). Heat is transferred through a groundwater flow system by conduction and convection. Conduct ive heat transfer occurs in the rock/water composite and is controlled by thermal conductivit ies of the porous medium and thermal gradients, while convective transport occurs under the influence of moving groundwater. Two types of convective transport are identified as convection and advection (also referred to in the literature as free and forced convection). The former term is related to the circulation of water due to forces caused by 3.3 Functional Relationships and Governing Equations 117 buoyancy, and the latter, the flow of heat and groundwater under a topographically driven system. The idea of using temperature data to estimate hydrogeologic parameters was first introduced by Bredehoeft and Papadopulos (1965). They showed how measurements of temperature wi th in a borehole could be used in a curve matching procedure to estimate vertical components of groundwater flow in the surrounding formation. M i n o r modifica-tions of the technique have been suggested by Stallman (1967) and Mansure and Reiter (1979). Examples of applying the technique to field situations are found in Sorey (1971) and Mansure and Reiter (1979). Examples of interpreting nonequilibrium temperature reversals measured in a series of boreholes to resolve a simple hydrogeologic structure and age of a geothermal system are given by Kasameyer et al. (1984), Blackwell (1985), and Eckstein et al . (1985). 3.3 Functional Relationships and Governing Equations Stal lman (1960) presented the governing equations for heat transfer in a saturated porous medium. More recently, Bear (1972) and Bear and Corapcioglu (1981) present developments of fluid, momentum and energy transfer in a thermo-elastic medium. The following equations are valid for a single-phase Newtonian fluid in a saturated porous me-dia. Solute concentrations are assumed to be negligible and cross-coupling phenomena such as the Soret and Dufour effect (Bear, 1972) are also assumed to be negligible. Normally, when considering deep regional flow systems, hydraulic conductivity shows significant departure from isothermal behaviour because of the dependence of fluid density and vis-cosity on temperature. However, for the temperature ranges considered in this study, these dependences are neglected. O n l y steady-state hydraulic head and temperature fields are considered in this study. This approximation is appropriate where yearly fluctuations in water table are small com-pared to the depth of the flow system and where seasonal temperature variations effect 3.3 Functional Relationships and Governing Equations 118 only the near-surface temperature field. Fluid Flow Thi s study wil l focus on the two-dimensional form of the steady-state groundwater flow equation: V • K • V * = uS{x - x') (3.1) subject to -K-VV -n = qf (3.2a) on Ti, and * = / ( * , y ) (3.26) on T 2 , V = (d/dx,d/dy), n is the unit outward normal, V - = (id/dx + jd/dy). Here, is the hydraulic head, OJ6(X - x') is a fluid source/sink term of strength u> at location x', K is the hydraulic conductivity tensor, qj is a specified fluid flux term on boundary T j , and f(x,y) is a function specifying Dirichlet boundary conditions on T2- If a numerical scheme is used to solve (3.1), a matrix equation results: h = A~1fi (3.3) where h is the approximate value of \l> due to the discretization, A is a global stiffness matr ix which is a function of K and grid design, and f\ is a loading vector wi th the appropriate boundary conditions. The fluid specific discharge is given by: q = - K • V * (3.4) Heat Transfer The appropriate heat transfer equation can be writ ten as: V • A - V 0 - pfCvfq-Ve = 0 (3.5) 3.3 Functional Relationships and Governing Equations 119 subject to - A • V 0 • n = qe (3.6a) on T3, and 0 = g[x,y) (3.66) on T 4 . Here, 0 is the temperature, pj is the fluid density, Cvf is the fluid specific heat capacity, A is the thermal conductivity tensor of the fluid/porous medium composite, q$ is a specified heat flux on boundary Ts, and g(x,y) is a function specifying Dirichlet boundary conditions on T4. If the external domain to T is also a porous medium, and fluid flow occurs across the boundary, then a component of thermal energy is carried wi th the moving groundwater. In this case, the appropriate boundary condition is of the third type. If the temperature on both sides of the boundary is assumed to be equal, then the boundary condition degenerates to V 0 • n = 0 (see Bear, 1972, p623). For the cases under consideration in this study, thermal conductivity is considered to be homogenous and isotropic; therefore A is treated as a scalar and (3.5) becomes: where D = p/CvfJ\ is a reciprocal-thermal diffusivity term. If a numerical scheme is used to solve (3.6), a matrix equation results: where 6 is the approximate value of 0 due to the discretization, B is a global stiffness matr ix which is a function of A , grid design and q, and J2 is a loading vector wi th the appropriate boundary conditions. W i t h i n the upper crust, heat is transferred primarily in the vertical direction. The equation for the vertical component of heat flow in a conductive environment is: V 2 0 - D g- V © = 0 (3.6) B = B~lh (3.7) (3.8) 3.3 Functional Relationships and Governing Equations 120 where qsz is the component of heat flow in the vertical direction. Variations in geothermal gradient [dO/dz] wi th depth are useful in interpreting and separating advective effects in thermal profiles. Non-linear Functionals In the inverse problem, the actual values of hydraulic heads and the model parameters are unknown. Usually (3.1) takes the form: h" = %x,p\ + v (3.9) where 3 is a non-linear functional relating k", the observed values of hydraulic head; x, a vector of grid coordinates in a numerical scheme; and p, the actual model parameters which could consist of hydraulic conductivities, boundary fluxes, and sources and sinks. v is a vector of residuals, which is assumed to have a mean of zero and some covariance, structure: E(u) = 0 E(vuT) = o2hVh (3.10) where o\ is a scaling factor used to determine the magnitude of the variances. If kriging is used to determine V^, then o\ is equal to one. Similarly, temperature is writ ten as: 6" = G{x,p)+$ . (3.11) where G is a non-linear functional relating the observed temperatures 0", and c is a vector of residuals assumed to have the following properties: £(<r) = o E(KT) = ojVe (3.12) where of is a scaling factor. If kriging is used to determine V#, then 0% is equal to one. Because the hydraulic head and temperature fields are measured wi th completely different 3.4 Construction 121 devices, it is reasonable to assume that there is no correlation between i> and f, so that E{ycT) = 0. However, if kriging were used to estimate both the hydraulic head and thermal fields, then there may be correlation between the two variables at interpolated points. A cross-covariance matrix may be developed in these cases. Output covariances of the estimated hydraulic head and temperatures from the inverse are also correlated. 3.4 Construction In Chapter 2, a construction philosophy was adopted for parameterizing a steady state flow model and a number of objective functions were presented to generate a variety of acceptable models. More confidence is placed on those features of the model that are replicated under all norms. In this paper general forms based on either the L 2 or L\ norms wi l l be used. The generalized L2 norm is: J = {h-h)TV-l(h- h:)-r \{Y - Y'fV-^Y - y * ) + 0{6- 6*)TV-l{0 - 0*) (3.13) where: , V y , and Vg are prior hydraulic head, log-conductivity, and temperature co-variance matrices, Y are the values of log conductivity determined by the inverse method, y * are observed or estimated log-conductivity values, and A and /? are scaling factors, which may be unknown. If 0", h*, and Y* are estimated by kriging, along wi th V^ , V y , and VQ, then A = 1 and 0 = 1. A or j3 may be treated as unknowns to be solved for in the opt imizat ion scheme if it is judged that the overall form of the prior covariance matrices are known to a greater degree than the magnitude of the errors. The inverse approach can be viewed as stochastic or deterministic depending on how the nature of the uncertainties in data, the model, and the noise vectors v or f are viewed (see Chapter 2, Appendix 2A) . Another group of objective functions under the L\ norm are suitable for construction purposes. Th i s norm is generally considered to be more robust than the L 2 norm as the underlying distributions in u and f could be viewed as being drawn from an exponential probabil i ty distribution function (pdf). Consequently, outliers in the data do not effect 3.5 I n v e r s i o n t e c h n i q u e 122 estimation of the model parameters to as great a degree as the L2 norm. These objective functions are: J = \h - h'\ + \\Y - Y'\ +(3\0 - 6'\ (3.14) and a weighted L\ norm: J = v ^ + ) r ^ i H / J r ^ ( 3 . 1 5 ) Here, M is the number of hydraulic head data, N the number of log-conductivity pa-rameters in the model, and O the number of thermal data. Also; o~hj, <7y, and o~gk are the standard errors of an observation in hydraulic head, log conductivity, and temperature, respectively. In the results section we wil l demonstrate the use of functionals (3.13) and (3.15). 3.5 Inversion technique The opt imizat ion technique described in Chapter 2 is used to solve the inverse problems posed in this study and the reader is referred to that chapter for a complete description of the method. In addition to computing the covariance structure of the estimated model cov pest, the covariances of the estimated hydraulic heads, cov hest, and the temperatures, cov 6est can also be computed. A t the termination of the optimization procedure, the temperature values at each node stored in a matrix T that correspond to each parameter set can be used to compute the covariances of the estimated temperature field. ! 1 {COV EEAT]L> = JL^T) ^ |T*FC ~^1|T^ " ^ (3-16) where > = i E T ^ ( 3 - 1 7 ) L 1 i t k = l and L is the number of search points in the simplex figure. Chapter 2 summarizes the stopping criteria for the optimization technique. In the con-strained simplex algorithm, when the average value of the objective function for al l search 3.6 Summary O f Algori thm 123 points is below the expected value of noise, the procedure is terminated and a sampling of the parameter hypersurface is carried out using all search points. The appropriate measure of misfit is the x 2 test if data errors are assumed Gaussian, independent with zero-mean. Here the x 2 misfit is defined as: X 2 = $ > , - - < ) 2 / " . ? (3-18) 1=1 In the constrained simplex method / is equal to M + O. If I > 5 the expected value of X 2 is approximately I and the standard error of x 2 is ox2 = \/2I. We define as acceptable any model with I - 0X2 < x2 < T + ox2 (3.19) The most likely model is the one corresponding to y 2 = I. Note that (3.13) can be writ ten in the form of (3.18) through the use of a spectral (eigenvalue) decomposition. Therefore, one may apply this statistical test to (3.18). In (3.18) we recognize that the expected value of the functional J is equal to o\[M -\- O), since the median of x 2 is approximately equal to the number of degrees of freedom, and the objective function is implici t ly scaled by a factor of o\. If the L\ norm functional is minimized (i.e equation 3.15) then a different strategy is adopted. This functional can be viewed as following exponential distribution. For the L\ norm the expected value of x 1 is ( f ) = + M) wi th a variance of (1 — ^) \0 + M] (Parker and M c N u t t , 1980). 3.6 Summary Of Algorithm In this section the computational strategy used in solving the inverse problem is sum-marized. 1. Decide on an appropriate interpretational model and parameterize the flow domain. Assign boundary conditions that are assumed to be known or treated as parameters to be solved for in the inverse. Discretize domain for the numerical solution to the 3.7 Inversion of Theoretical Data 124 boundary value problem. 2. Choose a functional that is to be minimized, for example: J — (h- h")TV~l(h - h*) -r\(Y - Y")TV~1(Y -Y') + 0(6 - 0*)TVe-\6 - $*) (3.19) Th i s functional can be written as: J = Jh + XJY + 0Je (3.20) Subject to: Yt<Y< Yu 3. Set the random number seed and randomly pick L trial points between the constraint boundaries. 4. Compute the hydraulic head at each nodal point. 5. Compute the specific discharges in each element. 6. Calculate the temperatures at each nodal point. 7. Compute functional J for each set of search points. The functional is computed by Cholesky decomposition as described in Chapter 1. 8. Follow the constrained simplex procedure unt i l the global min imum is reached at the appropriate noise level. 9. Compute the covariance of hest, pest, and 9est. 10. Check for non-uniqueness by resetting the random number seed and resolving the problem. 3.7 Inversion of Theoretical Data In order to investigate the robustness of the proposed approach, the technique is applied to a succession of two-dimensional problems. In some cases Gaussian random noise is added to theoretical data to simulate corruption of 'true' data values. Alternatively, errors in the 3.7 Inversion of Theoretical Data 125 data could contain structure, typical of kriging estimates. Only the former case is analyzed in the theoretical examples presented here. The following three statistics yield an estimate of the noise level in the data. These statistics are: 1 N M ^ • ^ j ^ E w - ^ W (3-21) i 1 se- = ( 7 E W - t f r T ) = (3-22) 3 = 1 The next set of three statistics quantify how well the inverse reproduces actual values of the data and model (haci , 6act , and Yact) in our theoretical examples: 1 N ^ = ( ^ E f t - y r T ) " (3-23) 1 M ^ = ( ^ £ f e - * n 2 ) * (3-24) ^—(yD^-^n 2)*. (3-25) where Y, h, 6 are the centroid of sample values of log-conductivity, hydraulic head, and temperature from the inverse algorithm. A s mentioned in the introduction to this chapter, there is a range to the advective disturbance for which this joint inversion scheme should show increased sensitivity of a change in parameter values to temperatures. Van der K a m p (1982, 1984) has shown how dimensionless numbers can be used to characterize advective effects. This author defines a B u l k Peclet number AW/H [ ' where Q is the total flow through the system per unit width, W is the total length of the flow field, and H is the total height of the flow field above a basal boundary. This 3.7 Inversion of Theoretical Data 126 number wi l l be presented along wi th model misfits when examining theoretical data. For values of P less than 1, the thermal regime is dominated by conductive effects. For values of P between 1 and 3, the thermal regime enters a transition and becomes advectively dominated. For P values much greater than 3, the thermal field is isothermal (Van der K a m p , 1982, 1984). Simple Two-dimensional Flow System In this first series of synthetic examples the ability of the joint inverse scheme to re-construct a simple structure is examined. The geometry of the flow domain, boundary conditions and log-conductivity distribution are shown on Figure 3.2. Two parameteriza-tions are examined. In Case A , a small inclusion of higher hydraulic conductivity material is placed in the middle of the flow system, and in Case B the size of the inclusion is ex-panded to cover a larger part of the flow system. The forward problem is solved by first superimposing a finite element mesh on the domain (Figure 3.2). The mesh consists of 77 nodes and 60 quadrilaterals (120 triangular elements). The finite element equations, (3.3) and (3.7), are solved wi th boundary con-ditions h = 1, 6 = 1 on one end and h = 0, 6 — 0 on the other. In a series of runs, the thermal disturbance induced by advective heat transfer is introduced by increasing the value of the reciprocal thermal diffusivity, D, from low (conductive) to high (advective) levels. This increase in D has an effect equivalent to modifying the Y\ and Y2 values by the same ratio of K\/K2 to induce higher specific discharge. The actual noise-free hydraulic heads and temperatures computed at each finite ele-ment node from the forward problem are used as the data set for the inverse algorithm. In these simulations it is assumed that the existence and location of the two zones is known w i t h certainty. /? is equal to one for all runs. Note that wi th hydraulic heads specified as boundary conditions, and if the thermal data are not included in the inverse, then the solution is non-unique. A n y Y\, Y2 values that preserve the correct ratio of K\jK2 w i l l generate identical hydraulic head values at all nodes. If thermal data are included 3.7 Inversion of Theoretical Data 127 in the inverse it is expected that within some range of fluid velocities, unique hydraulic conductivities can be resolved. In the heat transfer equation, temperatures are directly affected by the hydraulic conductivity of the medium. Note that only unique ratios can be determined in an inverse with Dirichet boundary conditions on the inflow and outflow boundaries. However, a simultaneous inverse scheme involving temperatures may resolve unique hydraulic conductivity values because the thermal equation is directly sensitive to absolute values in hydraulic conductivity. W i t h noise free data, the constrained simplex algorithm is terminated when five con-secutive iterations produce a difference between the maximum and minimum objective function values calculated for the search points that is less than a specified tolerance (Box, 1965). This strategy is adopted in noise free situations as computations of the data and model covariances are not applicable. Here, the opt imum value of the objective function should be close to zero, as the estimated and observed data can be closely matched. For Peclet numbers much less than 1, the thermal regime is dominated by the effects of heat conduction, so that the temperature field is not sensitive to the hydraulic parame-ters. In this case, the problem posed is non-unique. However, by starting the constrained simplex wi th the same seed number in each tr ial as the value of P is increased, the solu-tion converges to essentially the same set of parameter values for each t r ia l , yielding the consistent behaviour observed in the model misfit. The results of the joint inverse for Case A are shown on Figure 3.3 (top curve), for increasing P compared to the model misfit Sy. The range in P as shown indicates thermal fields ranging from conductive to advective-dominated levels. The top curve for Case A shows l i t t le change in model misfit unt i l l a value of P of about 1 is reached. A range of P between 1 and 10 is identified as the range at which log-conductivities show appreciable sensitivity to temperatures. M a x i m u m sensitivity is reached at about P equals 3. This sensitivity can be explained by examining the fluid and thermal fields. The temperature field for a conductive environment (Case A ) (P — 1 x 10~ 6 ) is shown on Figure 3.4a. A t 3.7 Inversion of Theoretical Data 128 higher fluid flow velocities (Figure 3.4b, P = 3.4), the 0.9 to 0.6 isotherms are clearly perturbed in the area of the inclusion, and are sensitive to the value of log-conductivity of this zone. A t sti l l higher flow velocities (Figure 3.4c, P = 11.2) all the isotherms are grouped at the right end of the flow system. The area containing the inclusion is essentially isothermal, and is therefore insensitive to the hydraulic conductivity of the inclusion. In contrast, the result of the inverse for Case B (Figure 3.3 bot tom curve) shows a more gradual improvement in model misfit for ranges of P up to about 3. This improvement in sensitivity can be attributed to the larger size of the inclusion. Here, larger portions of the thermal field are affected by the inclusion's size. As a result, the th i rd term of the objective function (3.13) involving temperatures has a greater effect on the minimizat ion. One further simulation is carried out with the small inclusion in Case A moved to the far right hand boundary in Figure 3.3 (not shown). Only a small decrease in SY is noted for all ranges of fluid flow. This apparent insensitivity of the inclusion on the right side can be at tr ibuted to the relative position of the inclusion wi th respect to the boundary conditions. Therefore, it can be concluded that the joint inverse scheme is sensitive to the Bulk Peclet number of the system, which is a measure of the level of advective disturbance. The inverse solution is also influenced by the position and scale of the hydraulic conductivity anomaly wi th in the flow field. Hypothetical Layered Medium In this next example the ability of the simultaneous inverse scheme to reconstruct a layered medium is examined. A common problem in groundwater hydrology involves determining the location and hydraulic properties of aquifers in a flow system. The example (Figure 3.5) is a hypothetical flow system of length 3000 m and a height of 1330 m. A low hydraulic conductivity layer at the base of the flow system from 0 to 300 m elevation restricts 'active' flow to the region wi th an elevation between 300 and 1330 m. B y imposing the basal layer wi th in the solution domain, any advective disturbance in the upper region 3.7 Inversion of Theoretical Data 129 is not constrained by an assumed conductive heat flux at the base of the active region of flow. Boundary conditions for the problem are also shown on Figure 3.5, and are chosen to be s imi lar to a later example from Downie Slide, Bri t ish Columbia. These conditions consist of a specified fluid flux of 1.2 x 1 0 ~ 9 m/sec applied along the recharge zone for the system, a hydrostatic head distribution on the far right hand side of 1320 m, and fluid flux from a point source located near the discharge end of the system (3.0 x 1 0 _ G m 2 / sec ) . The thermal boundary conditions are a specified temperature on the water table, and the left and right hand sides are adiabatic (a zero heat flux normal to the boundary). A heat flux of 108 m W / m 2 is applied at the base of the flow system. The thermal conductivity is 6.0 W / m ° C sec. Specifying a heat flux along this boundary instead of a specified temperature values results in a more realistic representation of the boundary condition. The hydraulic conduct ivi ty field consists of three values: a background value of 2.7 x 10~ 8 m/sec, a 20-110 m thick aquifer (from 870 to 980 m elevation at the far right hand boundary, Figure 3.5 ) with a hydraulic conductivity of 1.7 x 10~ c m/sec, and a basal layer with a hydraulic conduct ivi ty of 1.0 x 1 0 - 2 0 m/sec (log-conductivity values of -7.75, -5.75, and -20.00, respectively). Note the aquifer is not of uniform thickness. The problem is parameterized in terms of 11 layers of non-uniform thickness, each of which corresponds to a row in the finite element mesh. The sequence of steps following in this example begins by first solving the forward prob-lem using a finite element technique. The actual values of hydraulic head and temperature are corrupted by noise, and the inverse problem is solved using data from all nodal loca-tions. In subsequent simulations, information wil l be selectively removed and the ability of the inverse technique to reconstruct this hydrogeologic model wi l l be examined. The forward hydraulic and thermal problems are solved by first subdividing the prob-lem area into 310 rectangular elements (620 triangles) and then applying the finite element equations (3.3) and (3.7). Recall that both the groundwater flow system and the thermal 3.7 Inversion of Theoretical Data 130 regime are at steady state. Figure 3.6 shows the hydraulic head field and Figure 3.7, the resulting thermal field. The Bulk Peclet number calculated from the solution of the forward problem is 3.3. The inverse problem consists of determining the log-conductivity values in a number of layers in the system. If geologic layering is suspected prior to solving the inverse problem, one would likely parameterize the flow system in terms of layered zones. In this example, 11 layers are chosen. Each layer corresponds to a row in the finite element mesh used in solving the forward problem. The actual hydraulic head and temperature values at the 352 nodal points are corrupted with Gaussian random noise of N(0, 10 m) and N (0, 0.15 ° C ) , prior to solving the inverse problem. The upper and lower bounds on the log-conductivity parameters are specified as -5.00 and -9.00. No a priori data on hydraulic conductivity is included in the inverse, so the presence of the aquifer must be identified from information carried by the hydraulic head and temperature data. The lower hydraulic conductivity basal layer is assumed to be known and is fixed at -20.00. The L2 norm objective function (3.13) is minimized. Typical runs involve 6000 iterations and take about 1 hr of C P U time on a F P S 1 6 4 / M A X array processor. The results of this t r ia l are presented on Figure 3.8, which shows selected contours of log-conductivity , and in Table 3.1 (Case 1). This table shows significant improvements of the predicted hydraulic head and temperature data over the original noise corruption. Most significant is the small model misfit (Sy) value of 0.21 as a result of including the thermal data. Note on Figure 3.8 that the contour line for -7.75 wavers on either side of the aquifer. This meandering is a result of the noise corruption of the data and the plotting routines. Background values are determined to be close to a value of -7.75 and the joint inverse correctly identified the presence of the aquifer and determined a log-conductivity value in this layer of -5.76. In another t r ia l , the effects on the inverse of positioning the aquifer lower in the flow system is explored. In this case the aquifer is located between approximately 415 and 530 m elevation. The actual nodal values of hydraulic head and temperature are corrupted 3.7 Inversion of Theoretical Data 131 using the same noise level as in the previous example. A l l other parameters and objective function remain the same. The inverse correctly identified the more permeable zone wi th a log-conductivity value of -5.75 (Figure 3.9 and Table 3.1; Case 2). Background values in the other layers ranged from -7.5 to -7.9. On Table 3.1, the smaller value for the model misfit Sy suggests the solution is more sensitive for the deeper aquifer. The preceding simulations are based on a data set wi th a large number of values, each of which is coincident with a nodal location in the finite element mesh. In the next section the effects of a l imited data base are discussed, along wi th the effects of the uncertainties in boundary conditions. Hypothetical Layered Medium: Limited Data In the following example, an inverse problem wi th a more realistic data set is posed. A n inverse based on hydraulic head data alone wil l be carried out to identify the presence of an aquifer. This same exercise wil l be carried out again augmented with thermal data collected in the same observation wells as the hydraulic head data. The geometry of the problem and the location of an aquifer at depth are shown on Figure 3.10. The hydraulic conductivity field consists of three values: a background value of 2.7 x 1 0 - 8 m/sec, a 20-110 m thick aquifer (from about 415 to 530 m elevation) wi th a hydraulic conductivity of 1.7 x 10~ 6 m/sec, and a basal layer with a hydraulic conductivity of 1.0 x 1 0 - 2 0 m/sec (log-conductivity values of -7.75, -5.75, and -20.00, respectively). A fictitious exploration program is carried out to identify the hydrogeologic system. It is assumed that 16 boreholes are equally spaced across the flow domain. Borehole depths range from about 100 m in the discharge zone to about 500 m deep in the recharge zone. None of the boreholes penetrate the aquifer; hence the existence and hydraulic properties of the aquifer are not known. A n estimated value of hydraulic conductivity of 1.8 x 1 0 - 8 m/sec is determined for the porous medium. A tunnel is located in the domain, wi th a flow rate of 3.0 x 1 0 - 6 m 3 / s e c / m width. If the discharge remains reasonably constant over t ime, this suggests a steady-state groundwater flow system. The boundary conditions for 3.7 Inversion of Theoretical Data 132 the mathematical model of the flow system consist of the water table (position assumed to be known from dri l l ing), an approximate value of recharge to the system of 3.2 x 1 0 ~ 1 0 m/sec, and a hydrostatic head boundary at X equals 3 k m (based on regional groundwater flow considerations). There is a zone of low hydraulic conductivity at the base of the flow model . It is considered that this boundary is deep enough so that the mathematical solution of hydraulic head is sufficiently insensitive to variations in this parameter at the measurement points. It is further assumed that a total of 16 hydraulic head measurements are available across the system. These 16 values have noise levels of N |0 , 10 m], and correspond to completion zones at the bottom of the wells. The problem is parameterized in terms of 11 zones or layers as for the previous synthetic example. Each layer corresponds to a row in the finite element mesh. A n inverse based on the above data is carried out, minimizing a weighted L2 norm without prior information. Upper and lower bounds on log-conductivity are -9.0 and -5.0. Recharge rates to the system are also treated as an unknown to be determined in the inverse. Pr ior bounds on this parameter are 1 0 - 1 0 and 1 0 ~ 8 m/sec. The position of the water table is assumed to be known and is fixed in the inverse simulations. Note that an inverse based on hydraulic head head data set alone is potentially non-unique. Results of the inverse are tabulated in Table 3.2 and are shown in Figure 3.11. Table 3.2 (Case 1) shows the inverse determined a recharge rate of 2.8 x 1 0 ~ 9 m/sec, which is close to the correct value of 1.2 x 1 0 - 9 m/sec. The value of log-conductivity in the aquifer is determined to be -5.83 (true value -5.75). However, as shown on Figure 3.11, the log-conductivity values of the aquifer and layers immediately above and below are poorly resolved into one zone of Y > —6.75. This same exploration program is now augmented wi th temperature measurements from each of the boreholes. Temperatures can be profiled continuously down the holes, but only measurements that coincide wi th nodal coordinates in the finite element mesh are used. In this case 96 temperature measurements are available. It is assumed that the 3.7 Inversion of Theoretical Data 133 thermal conductivity of the domain is uniform and known. The basal heat flux is also assumed to be known and a ratio of 0.018 J C / m between these two quantities is fixed in the inverse. A simultaneous inverse of the combined thermal and hydraulic head data set is carried out using norm (3.13) without prior log-conductivity information. The results (Table 3.2: Case 2, Figure 3.12) show great improvement over the previous inverse using only hydraulic head data. A recharge rate of 1.2 x 10~ 9 m/sec is determined and the location and hydraulic conductivity of the aquifer are clearly determined. Sy values from the joint inversion are 0.22 , reduced from 0.64 from the hydraulic head inverse. Effects of Depth of Basin Smith and Chapman (1983) examine the thermal effects of restricting depth to low hydraulic conductivity basal layer. Perturbation of a conductive regime is greater for deeper flow systems. The main reason for this behaviour is that volume flow through the system increases for a basin with greater depth, assuming the permeability and water table posit ion remain fixed. As the volume flow increases so does the thermal disturbance (see equation 3.26). In some situations, when the location of a lower flow boundary is not well known a pr ior i , it may be necessary for the location of that boundary to be treated as another unknown parameter to be estimated. This determination amounts to solving for a hydraulic conductivi ty value of a lower basal confining layer, or set of layers, in the model. Let us return to the preceding data set. Previously a low hydraulic conductivity basal layer was fixed at 0 to 300 m elevation wi th a log-conductivity value fixed at -20.0. In one addit ional inverse simulation the hydraulic conductivity of this basal layer is determined. P r io r bounds on this parameter are set at -25.0 to -5.0. Results of the simulation show an identical result to the previous simulations (see Table 3.2, Case 2) and the log-conductivity value of the basal layer is estimated to be -17.0. Al though this value is inaccurate by 3 orders of magnitude, a low value of hydraulic conductivity is indicated, and thus the lower l imi t of the flow domain is determined. Addi t ional inverse simulations indicate that the value of hydraulic conductivity does not replicate well between simulations. 3.7 Inversion of Theoretical Data 134 This behaviour indicates a low sensitivity of variability of this parameter w i th respect to temperatures. Additional Synthetic Example In Chapter 2, an inverse analysis was carried out for a system where the hydraulic conductivities were generated from a stochastic model. Figure 3.13 shows this structure. The hydraulic conductivity values were generated from a lognormal distribution with a population mean of -7.00, and a standard deviation of 0.315. A n exponential covariance function was adopted, with an integral scale of 200 m . Recall that the log-conductivity field was first generated on a regular grid, and then the values were transferred to the finite element mesh. In this case the vertical correlation structure of the original population is not preserved. However, this simplification does not affect the purpose for which the example is designed. The steady-state flow problem is solved using the finite element technique. The problem area is sub-divided into 310 quadrilateral elements (620 triangles) with 352 nodes. The hydraulic heads at 33 points are sampled and corrupted with Gaussian noise of N(0,10 m). A residual kriging procedure is followed (see Appendix 2B) to interpolate the hydraulic head data to all nodes in the finite element mesh. A smallest deviatoric model can be constructed (equation 2.11) using hydraulic con-duct ivi ty values sampled from 31 points. In the inverse solution for this problem the water table position and all boundary conditions are assumed to be known. A n interesting com-parision can be made wi th these previous simulations and the joint inverse using thermal data. In these examples it is assumed that temperatures are recorded at discrete intervals wi th in a series of boreholes located at regular intervals across the basin (see Figure 2.11). The hydraulic head data set is augmented wi th 166 temperature data. The actual values of temperature are corrupted wi th noise of N(0,0.25°C) prior to the inverse solution. It is assumed that values of thermal conductivity and basal heat flux are known (6.0 W / m °C sec and 110 m W / m 2 , respectively). The Peclet number for this simulation is 3.0. Table 3.3 compares inverse results based on hydraulic head data alone (case 1), hy-3.7 Inversion of Theoretical Data 135 draulic head data with prior log-conductivity data (case 2), hydraulic head data with tem-peratures (case 3), and hydraulic head data wi th temperatures and prior log-conductivity data (case 4). a y is the standard deviation of the lognormal distribution computed from the inverse. Th i s number is compared to the ' true' value of 0.315. Examining the first column [Sf), one can see that the estimated hydraulic heads are fit to the about the same level of noise in the data in all cases. Comparing Case 2 to Case 3, the smallest devia-toric norm produces a better representation of the actual Y values than the joint inverse ( 5 y = 0 . 2 9 ) . However, the computed statistics of the realization from the joint inverse are better than the previous case (oy — 0.325). In other words, the joint inverse produces a structure wi th similar statistics as the actual field, but not in terms of closeness to actual values. Even wi th the inclusion of prior log-conductivity and thermal data (Case 4) the actual hydraulic conductivity field is no better determined than without thermal data. However, the benefits of including thermal data and applying a joint inversion can be substantial when considering the more realistic problem of uncertain boundary conditions. A simulation is carried out with an unspecified hydraulic head from C to C , an uncertain log-recharge rate from B to C , an uncertain basal heat flux, and an uncertain thermal conductivity. In this case the point sink of fluid is known as is the water table position and the temperature variation along the water table. This tr ial is labeled as Case 5 in Table 3.3. The computed values for oy of 0.33 indicate that the inverse replicates the actual standard deviation well. A value of Sy of 0.37 indicates that the parameter values obtained from the inverse with uncertain boundary conditions, are determined as well as if they are known with certainty. Estimates of the boundary condition parameters are given in Table 3.4. Good agreement between actual and computed values are indicated. Note that this same problem based on hydraulic head data alone or without a large base in prior hydraulic conductivity would be difficult to solve uniquely (see Chapter 2). It is our conclusion from the above simulations that the standard deviation of the original lognormal distribution is too low to perturb temperatures on a local (element by 3.8 Analysis of Field Data 136 element) scale. However, an inverse with temperature data does resolve the mean log-conductivi ty distr ibution. Note that changes in boundary conditions also affect the mean fluid-flow and therefore we are able to determine these parameters from a simultaneous inverse. 3.8 Analysis of Field Data In this section the inverse theory presented in this paper is applied a field problem. The site chosen is the Downie Slide in the east-central part of the Province of Br i t i sh Columbia, Canada . The slide has been studied extensively in the period from 1965 to present time. In Chapter 2, the regional hydrogeology, and water budget of the slide was reviewed. A hydraulic head data set was used in an inverse scheme to aid in the hydrogeologic interpretation of the slide. Results of Chapter 2 show that the only reasonable model resolvable at the site is that of a homogenous, isotropic medium. Four parameters were solved for; the free-surface position of the water table, and three boundary condtions for the domain. Only the ratio of the fluid flux to hydraulic conductivity could be resolved. It wi l l be demonstrated that wi th inclusion of thermal data, a value of homogeneous, isotropic hydraulic conductivity can be determined at this site. In the next few sections the thermal program at the slide is outlined, followed by the application of the simultaneous scheme to the combined thermal/hydraulic data set. Temperature Measurements Figure 3.14a shows a plan view of the slide with drillholes completed since 1965 and some of the principal instrumentation locations. The cross section that is modeled is shown on Figure 3.14b. Also shown are the locations at which temperature measurements have been recorded downhole. Temperatures in drillholes were recorded with a thermistor probe connected to a portable Hewlett Packard 3465A 5-place digital voltmeter. Y S I thermilinear thermistors 3.8 Analysis of Field Data 137 were used in the thermal probe. These thermistors have a nominal resistance of 10,000 ohms at 20°C, operating temperatures of -50 to +50 °C , and a response time of 1 s. The probes were calibrated in a thermal 'bath ' at the University of Br i t i sh Columbia. The kerosene bath was used to compare the probe thermistor and resistance to that of an Omega pla t inum standard thermometer, by determining probe resistance R and tem-perature O over a range of 1 to 17 ° C . Temperatures of the standard were recorded by measuring the standard's resistance wi th a Fluke 88O0A 6-place digital voltmeter and con-verting to temperature wi th calibration charts. This procedure produces a series of probe-resistance (R) vs temperature points for the field probe. Rather than fit a line or smooth function through the calibration points, field probe resistance readings were converted to temperature by linear interpolation between calibration points. The accuracy of a field temperature measurement is limited by accumulation of calibration errors and stability of the resistance meter, and is estimated to be ± 0.1 I J C . In the field study, measurements were made at discetete depth intervals (commonly 5-10 m ), after wait ing at least one minute for the probe temperature to equilibrate wi th its surroundings. Measurements were taken both above and below the water table, however, only those measurements taken below the water table are used in the inverse. The field program was carried out in July-August , 1985, and temperature-depth profiles for 13 sites are shown on Figure 3.15. The drillholes vary in depth from 133 to 259 m. The drillholes are grouped into 4 zones based on collar elevation on the slide: Up-slope (S8), Mid-slope (S3, S12, S31, S34, S36), Bottom-slope (S14, S30, S37, S U , S32, S40), and Toe (1-04-01C, S17). Thermal profiles in holes located in recharge zones typically exhibit convex warping due to advective movement of cooler waters to deeper regions (Smith and C h a p m a n , 1983). The Up-slope, M i d , and Bottom-slope temperature profiles show this behaviour. Based on the hydrogeologic modeling discussed in Chapter 2, this behaviour is reasonable as downward flow is prevalent in the near-surface hydrogeologic regime. 3.8 Analysis of Field Data 138 Thermal Conductivity Determination Unlike hydraulic conductivity which can vary over 13 orders of magnitude, thermal conductivi ty in rocks varies by less than a factor of 30. (Clark, 1966). The primary factors which cause variation in most silicate rocks are differences in quartz content, porosity and temperature. A useful equation that duplicates laboratory results for relating thermal conduct ivi ty of a rock matr ix to that of a porous medium is: A = a £ - * a £ (3.27) (Sass et al . , 1971); where A& is the thermal conductivity of the rock matr ix, A„, is the thermal conductivity of water (0.586 W / m sec °C at 10 0 C) and 4> is the porosity. On Downie Slide poTosity is estimated to be 1 % (Freeze, 1982). Thus the second term in the above equation is not likely to significantly affect the bulk A term. O n Downie Slide the thermal conductivity of the rock mass is considered homogeneous for the purposes of mathematical modeling. This assumption is based on the uniformity of the geologic materials at the site and the lack of coherent stratigraphy on the slide (see Chapter 2). Ranges in thermal conductivity for metamorphic rocks are between 1.5 to 4.7 W / m sec ° C , w i th a mean value of 3.3 (Clark, 1966). A best fitting thermal conductivity value wi l l be determined from the inverse. Thermal measurements in boreholes S32 and S30 provide an indication of the magni-tude of conductive disturbances due to thermal conductivity variations. These two holes have inclinometer casing fully grouted to the borewalls over their entire length (Figure 3.15). Large scale variations such as convex or concave warping are due to advective ef-fects from the groundwater flow system. Small scale variations (less than ± 0 . 2 ° C ) are likely due to local variations in thermal conductivity or porosity. Thermal Disturbances M a n y factors affect subsurface temperatures and geothermal gradients measured at a given location. The types of effects which can be important in temperature measurements 3.8 Analysis of Field Data 139 for individual holes are: 1 Conductive effects associated with temperature or thermal property variations along the surface or at depth, 2 Advect ive disturbances due to subsurface water flow, either induced by the drill ing or existing on a regional or local basis. The dri l l ing process can cause significant disturbance of the temperature in the rock around the hole because of circulation of cool dri l l ing muds. Because of the large difference in time between the dri l l ing of the holes and the temperature measurements, these effects can be ignored. Transient surface-temperature variations are known to penetrate to depth in conductive environments. Diurnal surface variations (maximum perturbation of A © ± 0.01°C) are about 0.5 to 1 m for daily variations and about 20 m for annual variations (Gretner, 1981). Glaciat ion and long-term climatic effects penetrate to much greater depths (500 m). In some Pleistocene and Holocene sediments these variations can lead to corrections of 5 to 20 % higher surface heat flow than unadjusted values (Gretner, 1981). For advective environments such as Downie Slide, these short and long term climatic effects can be ignored (Lachenbruch and Sass, 1977). Short term effects can neglected due to the depth to water table below surface while long term effects are not considered because the advective transfer of heat damps such effects over a relatively short time frame. Borehole (casing) effects can be grouped into internal and external effects. Internal heat effects can be caused by thermal convection due to bouyancy of the fluid within the borehole itself. Holes over 7.6 cm may show convective effects. These effects are usually el iminated by installing smaller (2.5 - 5 cm) diameter casing wi th in the borehole. The observations at Dowine Slide are carried out in 1.9 cm P V C pipe or 5.9 cm I.D. A B S plastic or 7.5 cm I .D. aluminum inclinometer casing, so the readings are not considered to 3.8 Analysis of Field Data 140 be affected by these phenomena. External casing effects are caused by downward or upward flow of water within the hole or in the annulus between the casing and the borehole wall . This phenomenon has been recognized for some time and is a major reason that most holes used for thermal studies are cased or grouted following completion. If higher or lower hydraulic conductivity zones exist in the rock mass, hydraulic heads may be substantially different in each of these zones. Al though flow wi th in the rock mass may be such that different zones are not connected by fluid flow with a high enough rate to perturb the conductive regime, the borehole 'short-circui ts ' less permeable zones and allows for rapid intra-zone flow. As a consequence, water wi l l flow up or down sections of the hole affecting temperatures wi th in the flow zone and in some cases above or below the entry points of flow. This effect can be recognized on temperature/gradient depth plots. Usually, an isothermal zone on a plot has high or low gradients at either end of the zone. O n Downie Slide, boreholes S17 and S8 are affected over various depths. In S8, all but the bottom-most portion of the hole is affected by what is probably downward flowing cooler water in the annulus between the inclinometer casing and the wall of the hole. This hole was one of the earlier installations (1965) and the annulus is believed to be only partially filled with sand. Other sections may be void of backfill materials or contain clay-gouge material at shear zone locations. Downward flow is expected wi th in the general vicini ty of S8, because the area is a recharge zone for the groundwater flow system (see Chapter 2). S17 is one of the newer (1974) inclinometer holes, but it was only part ial ly backfilled with fine sand. Intra-zone water flows were encountered here during dri l l ing, so that upward moving groundwater is expected at this locale. Only the top 90 m is considered to be affected by rising warmer groundwater. The bot tom portion from 90 to 183 m is considered free of casing effects. The remainder of the holes at Downie Slide are either fully-grouted inclinometers or grouted (in zones) multi-level standpipe piezometers. Most of the multi-level standpipes are considered to be completed such that piezometer completion zones are adequately blocked from one another, 3.8 Analysis of Field Data 141 thus reducing or eliminating casing effects. This conclusion is supported by examination of gradient depth plots. Assumptions in the Analysis Chapter 2 describes the assumptions used in applying an inverse model to this site using hydraulic head data. The assumptions include: two-dimensional steady flow, and the ho-mogeneous, isotropic nature of hydraulic conductivity . In the joint inverse scheme, these assumptions are made, and also, additional assumptions regarding the thermal boundary value problem must be made. These assumptions are: 1. Two-dimensional approximation to the thermal field. Here the surface relief in a third direction is small compared to the relief of the two-dimensional section. The topogra-phy of Downie Slide and surrounding environs is typical of this situation. Downie Slide is a geologic feature on the old Columbia River Valley. This larger feature forms an extensive trough which extends for hundreds of kilometers. Smaller tributary valleys empty into the Revelstoke Reservoir near the slide from both sides, but on the large scale, the topography of the valley could be catagorized in terms of two dimensional relief. In addition, the effect of a fluid sink (drainage adit system) on the thermal field is to hold the water table position in a relatively constant position in the third direction. 2. Steady-state thermal field. The thermal field is subject to transient responses due to the seasonal nature of air temperature. The data used in the inverse is collected well below the depth of these transients. For an advective environment such as Downie Slide, a steady-state thermal field is reasonable. 3. Projection of boreholes off-section. Because of projections, measurement points wi l l not be exactly at the same elevation on different sections. This wi l l introduce an unknown component of noise into the measurements. These sources of noise are probably small in comparsion to other assumptions. 4. Uncertain boundary conditions. On the section under consideration there are uncer-3.8 Analysis of Field Data 142 tainties in hydraulic head, temperature and flux quantities on most of the boundaries of the domain. The various boundary conditions, either assumed or unknown, are shown on Figure 3.16. The left side of the domain is assigned as an insulated boundary due to symmetry on the adjacent side of the valley. Here fluid flow is forced to be vertical and this condition leads to vertical heat flow. A similar condition is applied to the far right boundary at X equals 3 km (C - C ) . Here a specified hydraulic head boundary (near hydrostatic) is assumed (see Chapter 2) and across this boundary moving fluid transports thermal energy. The appropriate thermal condition along this interface is zero heat flux normal to the boundary (see Bear, 1972, p 623). Along the water ta-ble a specified temperature boundary allows for both conductive and advective heat transfer. One would ideally like to know this condition along the entire water table. In mountainous terrain a reasonable lapse-rate of temperature change wi th elevation is 7 ° C / k m . Unfortunately, there is not enough thermal data to accurately assign wa-ter table temperatures. Furthermore, the water table position is unknown a priori . Therefore, the temperature along the water table is left as an unknown to be solved for in the inverse. The drainage system at the slide is approximated as a point sink of fluid which also carries an advective quantity of heat out of the system. The equivalent thermal condition for this point is left as a parameter for which the temperature is determined from the inverse. Along the base of the model an unknown conductive heat flux is applied across an assumed low hydraulic conductivity layer of 300 m thickness. The value of heat flux along this surface is assumed to vary linearly between both ends. In this way allowance is made for the possibility of uneven heat flow. The assumption of the known hydraulic conductivity value of the basal layer wi l l be relaxed in some simulations. The value of the basal heat flux is left as an unknown to be determined in the inverse. Basal heat flux values should be in the range of 100 to 120 m W / m 2 , based on heat flow studies south of the project area (Lewis, 1985, personal communication). 3.8 Analysis of Field Data 143 Application of The Inverse Technique A n ini t ial step in applying the inverse procedure in aquifer problems is often to use kriging or some other method to interpolate the values of observed data to all nodal points in a finite-element system (for example, Clifton and Neuman, 1982). The difficulties in apply ing interpolation techniques to the hydraulic head data at this site were examined in Chapter 2. The problems of interpolating thermal data are the same as for hydraulic head data; namely that a drift component exists in the data, and most of the data points are collected near the water table and are concentrated vertically. These factors make the application of kriging difficult. Interpolation of the temperature field was attempted by a residual kriging technique (see Appendix 2-B) . To test whether or not the resulting temperature field is unbiased, trials have been carried out using the thermal field generated from the last synthetic example. A vertical orientation of measurement points is assumed, s imilar to the work in Chapter 2. Based on the results of the kriging exercise, it is not possible to generate an unbiased estimate of the thermal field to al l nodal points in the finite element mesh. Thus, no attempt is made to use kriging on the temperature data set collected at the slide. There are limitations in applying the joint inverse method to Downie Slide. Whi le a total of 13 boreholes have been thermally profiled during the field program, only 7 of these are near to the section B - B ' . The requirement that a finite element node must be coincident w i th a thermal data point places an upper l imit on the number of data points available to the inverse because of computer storage limitations. Also, some boreholes on section have portions of the borehole thermally disturbed by casing effects. There are also uncertainties in boundary conditions and domain geometries in the form of a free-surface over about one-third of the water table surface. T h e simplest model that can be resolved from the data base is a homogeneous, isotropic hydraulic conductivity structure, wi th boundary conditions and thermal properties. This inverse can be classified as deterministic (constrained-weighted least squares). 3.8 Analysis of Field Data 144 The various parameters are shown on Figure 3.16 with specific prior ranges detailed in Table 3.5. The parameters to be solved for in the inverse are as follows: 1. The free-surface position of the water table; A to C represents the free-surface location, with the input flux parameter applied from B to C 2. The log-conductivity value of the basal layer, 3. The log-conductivity values of the landslide mass, 4. The specified hydraulic head value at the base of the flow system at C A n adjustable constant head boundary is applied along C - C . 5. The log-flux recharge rate to the system from B to C 6. The adit sink-strength (m 3 / sec /m) 7. Basal heat flux ( m W / m 2 ) 8. The thermal conductivity of the domain ( W / m ° C sec) 9. The temperature at the adit. In addition, the temperature of the water table as a function of position is another parameter (not shown). The basic procedure as noted under "Summary of Algor i thm" is followed with the exception that forms of (3.13) and (3.15) without log-conductivity penalty terms are min-imized. A n inner-iteration loop is defined, in which the free-surface is adjusted under each parameter set unti l the elevation of the surface equals the hydraulic potential. The matri-ces Vh and Vg are taken as a diagonal matrices. A s noted by Neuman and Yakowitz (1979), when the number of observation points is l imited an accurate estimate of the covariance mat r ix Vh or V$ may become difficult. In such cases o\ and a 2 may be unspecified, and this corresponds to a situation where Vh and VQ represent relative weights between ob-servations, and the overall scale of a priori error is unknown. Experience wi th the joint inverse shows that the choice of an a priori o\ or a\ does not influence the centroid of the parameter estimates to a large degree; here the o2 terms only effect the scale of covariances of hest and 0est . The results of the joint inversion runs are shown on Table 3.6. Typ ica l runs took 3.8 Analysis of Field Data 145 about one hour on a F P S 1 6 4 / M A X array processor. Results presented in Table 3.6 are averages of the parameter sets determined from the inverse. Each of the t r ia l sets produced data which satisfies the observations. Similar results between two norms indicate a strong possibil i ty that the final parameters are good representors of the 'true' model. Several runs w i t h different norms and different inital simplex configurations are presented; of interest is the spread in model parameters for some parameters under different norms and different starting configurations, indicating that some of the model parameters are poorly resolved. For instance, the basal heat flux, the thermal conductivity, and the specified head at C vary by relatively large amounts between simulations. These variations are due to an insufficient data base to adequately resolve these parameters. None the less the constructed models do offer valuable insight into the hydrogeology of the field site, which is not apparent from the inverse based on the hydraulic head data set alone. The constructed models persistently show a hydraulic conductivity structure of about 1 x 10~ 7 m/sec, which is consistent with previous estimates of hydraulic conductivity at the site (see Chapter 2). The basal layer log-conductivity is determined to be -8.6, which does indicate that a low hydraulic conductivity base in the flow system is adequately determined. A further comparison with the inverse results in Chapter 2 show good agreement between the two inverses for the hydraulic properties. Comparing ratios of recharge rates to hydraulic conduct ivi ty (qj/K) and adit sink strength to hydraulic conductivity ( w / K ) with previous results (Table 2.9) show very similar quantities. Recharge rates are estimated to be about 7 x 1 0 - 9 m/sec; close to the estimated value of 1 x 1 0 ~ 8 m/sec based on meteorological studies and baseflow calculations (see Chapter 2). Basal heat flux values are found to be in the range of 100 to 114 m W / m 2 , which is reasonable based on heat flow studies south of the project area (Lewis, 1985, personal communication). Thermal conductivity values are determined to be about 5.0 W / m °C sec; a value that indicates a high silica content in the rocks, more typical of quartzitic gneisses than schists. The temperature field determined from the inverse ( $ e s t ) is shown on Figure 3.17, wi th water table temperatures determined 3.9 Discussion and Conclusions 146 along the free-surface. A warping of temperatures towards the adit is apparently due to increased velocities at this area. 3.9 Discussion and Conclusions In the solution of the inverse problem, the ability of the data to resolve the model may be poor. One way of overcoming data limitations in an inverse problem is to intro-duce a second data set for another potential and apply simultaneous inversion. The joint inversion method seeks to use data from a number of different measurements to improve the resolution of parameters which are in common to one or more functional relationships. One such data set is subsurface temperature, which is sensitive to variations in hydraulic conduct ivi ty in more permeable systems. The nature of an advective disturbance has a significant impact on the proposed method. In low-permeability systems, where the thermal regime is conductive, inver-sion of the temperature field to construct hydraulic conductivity models wi l l not yield good results; the temperature field being far too insensitive to hydraulic conductivity vari-ations. Similarly, at high fluid velocities variations in hydraulic conductivity wi l l not effect temperature patterns. The basic idea behind the joint inversion scheme is simple; temperature measurements are taken in piezometers or boreholes along wi th hydraulic head measurements. Many more temperatures than hydraulic head measurements can generally be taken in most situations. This feature increases the amount of temperature data compared to hydraulic heads that can be gathered at most sites. In addition, temperature measuring equipment is s imple to construct and use in the field. The simultaneous inverse scheme is tested on a variety of synthetic examples. It can be concluded from studies of simple hydraulic conductivity structures that the joint inverse scheme is sensitive to the Bulk Peclet number of the system, which is a measure of the level of advective disturbance. In other examples the ability of the simultaneous inverse scheme 3.9 Discussion and Conclusions 147 to reconstruct a hypothetical layered medium is examined. One of the common problems in groundwater hydrology is determining the location and hydraulic properties of aquifers in a flow system. This problem is particularly difficult if exploration boreholes do riot penetrate the aquifer or where boundary conditions are uncertain. W i t h a limited data base in hydraulic head the inverse indicates a smearing of the hydraulic conductivity in zones below measurement points. In the synthetic examples chosen, the inclusion of thermal data correctly identified the recharge rate and the location and hydraulic conductivity of the aquifer. In another series of simulations the lower boundary to a flow system is determined wi th the joint inverse. Addi t iona l runs on a heterogenous medium presented in Chapter 2 have been carried out. W i t h a good temperature data base, thermal properties can be properly resolved. However, this stochastic problem (with certain boundary conditions) the addition of ther-mal data did not condition the inverse to a greater degree than accomplished with the addition of prior information on log-conductivity. The benefits of including thermal data and applying a joint inversion can be substantial when considering the more realistic problem of uncertain boundary conditions. A simulation is carried out with unspecified hydraulic head boundaries, an uncertain recharge rate, uncertain basal heat flux and an uncertain thermal conductivity. Good agreement between actual and computed values is indicated. Th i s same problem based on hydraulic head data alone or with a large base in prior hydraulic conductivity data would be difficult to solve uniquely. The simultaneous inverse is also applied to the Downie Slide data set examined in Chapter 2. Results from that chapter show a simple model of a homogenous, isotropic medium. Four parameters are solved for: the free-surface position of the water table, and three boundary conditions for the domain. Unfortunately, wi th a homogeneous hydraulic conductivity, al l that can be determined from the inverse are ratios of flux to hydraulic conductivity. B y including thermal data, the value of hydraulic conductivity can be deter-mined at this site. Some of the model parameters (basal heat flux, thermal conductivity, 3.9 D i s c u s s i o n a n d C o n c l u s i o n s 148 specified head boundaries) are not resolved well by the joint scheme. None the less the constructed models do offer valuable insight into the hydrogeology of the field site. The constructed models persistently show a hydraulic conductivity value of about 1 x 10~ 7 m/sec, which is consistent wi th previous estimates of hydraulic conductivity at the site. A basal layer log-conductivity is determined to be -8.6, which indicates that a low hy-draulic conductivi ty base in the flow system is inferred. A further comparison wi th the inverse results in Chapter 2 show good agreement between the two inverses for the hy-draulic properties. Recharge rates are estimated to be about 7 x 1 0 ~ 9 m/sec; very close to the estimated value of 1 x 10~ 8 m/sec based on meteorlogical studies and baseflow calcu-lations. Basal heat flux values are found to be in the range of 114 to 100 m W / m 2 , which is reasonable based on heat flow studies south of the project area. Thermal conductivity values are determined to be about 5.0 W / m ° C sec: a value that indicates a high sil ica content in the rocks, more typical of quartzitic gneisses than schists. Notation A global stiffness matrix for fluid flow B global stiffness matrix for heat flow Cvj fluid specific heat capacity cov he3t estimated hydraulic head covariance matr ix cov h' a priori hydraluic head covariance matrix cov peat estimated model covariance from inverse cov Best estimated temperatures from inverse functional transformation / function specifying Dirichlet boundary conditions (hydraulic) / i load vector for fluid flow equation /2 load vector for heat transfer equation G functional transformation function specifying Dirichlet conditions (temperature) H total height of the flow system Storage matr ix : k ' th estimate of i ' th value of hyd: raulic head h hydraulic heads computed by finite element. h* observed or interpolated (from kriging) hydraulic heads hest estimated hydraulic heads from inverse, based on mest h centroid of hydraulic heads from inverse K average value of the i ' th hydraulic head I identity matr ix / number of data points J objective function K hydraulic conductivity tensor L lower triangle matrix in Cholesky scheme L number of search points Chapter 3 - Notation 150 • L\, Li vector norms referring to functional types M number of observed and estimated data points m N dimensional vector of parameters N number of parameters O number of temperatures q specific discharge of fluid qj specified fluid flux on the boundary Sh root mean square deviations from inverse to actual data Sh~ root mean square deviations from observed to actual data Sy root mean square deviations from inverse to actual value Sy root mean square deviations from observed to actual value T*k Storage matr ix: k ' th estimate of i ' th value of temperature V scaled covariance matrix of residual heads Vh scaled covariance matrix of estimated hydraulic heads VY scaled covariance matr ix of estimated log-conductivity Ve scaled covariance matrix of estimated temperatures W total length of flow system x vector of grid coordinates in a numerical scheme Y, Y, Y' log-hydraulic conductivity, centroid of all Y values from the inverse, and observed or interpolated (from kriging) log-conductivity Yi, Yu prior constraints on log-conductivity in the form of lower and upper bounds 6 hydraulic heads computed by finite element. 0' observed or interpolated (from kriging) hydraulic heads 0e3t estimated hydraulic heads from inverse, based on mest e centroid of hydraulic heads from inverse average value of the Tth hydraulic head r i , r 2 boundary domain segments Chapter 3 - Notation 151 7 integral scale of a random variable A scaling parameter in multi-objective criteria 0 scaling parameter in multi-objective criteria £ vector of residuals between computed and observed temperatures A thermal conductivity of the porous medium v vector of residuals between computed and observed hydraulic head theoretical hydraulic head Pj density of the fluid u sink strength of adit Oh scaling factor used in defining magnitudes of covariances in hydraulic head OQ scaling factor used in defining magnitudes of covariances in temperatures oy scaling factor used in defining magnitudes of covariances in log-conductivity o~xi s tandard deviation of the chi-squared variable X 2 chi-squared statistic X 1 chi-one statistic C h a p t e r 3 - References 152 References Bear, J . , 1972. Dynamics Of Fluids In Porous Media. Elsevier, New York. Bear, J . and M . Y . Corapcioglu, 1981. A mathematical model for consolidation in a ther-moelastic aquifer due to hot water injection or pumping, Water Resources Research, 17(3 ) , 723-736. Blackwel l , D . D . , 1985. A transient model of the geothermal system of Long Valley Caldera, Cal i fornia , J. Geophys. Res., 9 0 ( B 1 3 ) , 11229-11242. Bodvarsson, G . , 1973. Temperature inversions in geothermal systems, Geoexploration, 11, 141-149. Box , M . J . , 1965. A new method of constrained optimization and a comparison with other methods, Computer Journal, 8, 42-52. Bredehoeft, I D . and l .S . Papadopulos, 1965. Rates of vertical groundwater movement estimated from the earth's thermal profile, Water Resources Research, 1(2), 325-328. C la rk , S.P., J r . , 1966. Thermal conductivity. In Handbook Of Physical Constants, (S.P. Cla rk . Jr . , ed.), Geol . Soc. A m . M e m . 97, 459-482. Cl i f ton , P . M . and S.P. Neuman, 1982. Effects of kriging and inverse modeling on condi-t ional simulations of the A v r a Valley aquifer in Southern Arizona, Water Resources Research, 22(2), 1215-1234. Eckstein, Y - , G . Maura th and R . A . Ferry, 1985. Model ing the thermal evolution of an active geothermal system, J. Geodynamics, 4, 149-163. Freeze, R . A . , 1982. Groundwater conditions on Downie Slide, 1, Review of the available hydrogeological data on Downie Slide. Report to Hydroelectric Design Division, Bri t ish Co lumbia Hydro , Vancouver, Canada. Gretner, P . E . , 1981. Geothermics: Using temperatures in hydrocarbon exploration. Edu-cation Course Notes Series Meeting, San Francisco. Kasameyer, P . W . L . W . Younker and J . Hanson, 1984. Development and application of a hydrothermal model for the Salton Sea Geothermal F ie ld , California, Bull. Geol. Soc. Amer., 9 5 ( 1 0 ) , 1242-1252. Lachenbruch, A . H . and J . H . Sass, 1977. Heat flow in the United States and the thermal regime of the crust. In The Earth's Crust, Geophys. Mon. Ser., 20, ( J .G . Heacock, ed.), A m . Geophys. Un ion , Washington, D . C . . Mansure, A . J . and M . Reiter, 1979. A vertical groundwater correction for heat flow, J. Geophys. Res., 8 4 ( B 7 ) , 3490-3496. Neuman, S.P. and S. Yakowitz , 1979. A statistical approach to the inverse problem of aquifer hydrology, 1, Theory, Water Resources Research, 15(4) , 845-860. Chapter 3 - RofVicurrs 153 Parker, R . L . and M . K . M c N u t t , 1980. Statistics for the one-norm misfit measure. ./. Geophys. Res., 8 5 ( B 8 ) , 4429-4430. Reklai t is , G . V . , A . Ravindran and K . M . Ragsdell, 1983. Engineering Optimization. Meth-ods and Applications. Wiley, New York. Sass, J . H . . A . H . Lachenbruch and R . J . Monroe, 1971. Thermal conductivity of rocks from fragments and its application to heat flow determinations, J. Geophys. Res.. 7G. 3391-3401. Smi th , L . and D.S. Chapman, 1983. O n the thermal effects of groundwater flow, 1 Regional scale systems, J. Geophys. Res.. 8 8 ( B 1 ) , 593-608. Sorey, M . L . , 1971. Measurement of vertical groundwater velocity from temperature profiles in wells, Water Resources Research, 7(4) , 963-970. Sta l lman, R . W . , 1960. Notes on the use of temperature data for computing groundwater velocity. In 6th Assembly on Hydraulics, (Soc. Hydrotech. France), Nancy, June 28-30. Stal lman, R . W . , 1967. Flow in the zone of aeration. In Advances in Hydroscience, V. 4 ? (Ven Te Chow, ed.), Academic Press. van der K a m p , G . , 1982. Interactions between heat flow and groundxvater flow - A review. Pro j . 109-17, Waterloo Res. Int., Un iv . of Waterloo, Waterloo. Ont.. van der K a m p , G . , 1984. Evaluating the influence of groundwater flow systems on geother-mal conditions. In Proc. Energ. '84, (Regina). Saskatchewan. Woodbury, A . D . and L . Smith, 1985. O n the thermal effects of three dimensional ground-water flow, J. Geophys. Res., 9 0 ( B 1 ) , 759-767. Ziagos, J .P . and D . D . Blackwell , 1981. A model for the effect of horizontal fluid flow in a thin aquifer on temperature-depth profiles, Trans. Geotherm. Resour. Counc, 5, 221-223. Chapter 3 - Tables 154 Table 3.1 Synthetic hydraulic head and temperature data set. Sh- = 10.40, S§ = 0.16. L2 norm results. Case 1: aquifer at mid-basin level, Case 2: aquifer at depth. Case Sh Sfi SY 1 1.73 0.03 0.21 2 2.70 0.02 0.07 Chapter 3 - Tables 155 Table 3.2 Sparse synthetic hydraulic head and temperature data set. Statistics of log conductivity dis t r ibut ion: S V = 10.40, S^= 0.16. Value of log-flux = -8.92, value of Y in aquifer = -5.75. Li norm results. Case 1: 16 hydraulic head data points only, Case 2: data set agumented wi th 96 temperature measurements. Case Sh S§ SY log — flux Y aquifer 1 13.38 - 0.64 -8.55 -5.83 2 9.77 0.11 0.22 -8.91 -5.74 Chapter 3 - Tables 156 Table 3.3 Synthetic hydraulic head and temperature data set. Statistics of log conductivity realiza-t ion: HY = —7.00, OY = 0.315, S^- — 7.63, Sg- = .460. L2 norm results. Case 1: hydraulic head data only, no prior information, Case 2: hydraulic head data with 31 prior Y data, smallest deviatoric norm, Case 3: hydraulic head data wi th 166 temperatures, Case 4: hydraulic head data wi th 166 temperatures and 31 prior Y data, smallest deviatoric norm, Case 5: hydraulic head data wi th 166 temperatures and uncertain boundary conditions. Case Sh S§ E(Y) o~Y Sy 1 10.1 - -6.95 .41 .48 2 9.9 - -6.98 .18 0.29 3 11.5 .70 -6.98 .33 0.42 4 12.4 .88 -6.99 .16 0.30 5 10.9 .74 -6.97 .33 0.37 Chapter 3 - Tables 157 Table 3.4 Comparison of actual and inverse results. Case 5: hydraulic head with 166 temperatures and uncertain boundary conditions. Statistics of log-conductivity distribution: ny — - 7 . 0 0 , ay = 0.315, Sh- = 7.63, Se- = .460. L2 norm results. Description BC1 BC2 BC3 BC4 Actua l 1333 -8.01 108 6.0 Inverse 1333 -8.14 104 6.2 B C 1 Value of specified hydraulic head at X equals 3000 m , Y equals 0 m. B C 2 Log-recharge rate (input fluid flux) B C 3 Basal heat flow B C 4 Thermal conductivity Chapter 3 - Tables 158 Table 3.5 Input parameters and ranges: Downie Slide joint inversion Parameter Lower Bound Upper Bound Y - Basal Layer -20 -5 Y - Slide Mass -9 -5 Specified Head (m) C 1333 1450 Log-recharge rate (flux) -10 -7 Log-sink strength -4.7 -4.5 Basal Heat flux ( m W / m 2 ) 90 120 Thermal conductivity ( W / m °C sec) 1.5 6.5 A d i t temperature °C 6 14 Chapter 3 - Tables 159 Table 3.6 Downie Slide inverse results: joint inversion Parameter L\ Y - Basal Layer -8.17 -8.21 -8.13 -8.20 Y - Slide Mass -7.11 -7.08 -7.02 -7.06 Specified Head (m) C 1375 1367 1380 1347 Recharge rate (flux) -8.17 -8.24 -8.09 -8.05 Log-sink strength -4.62 U. 62 -4.57 -4T61 105 - 10H Basal heat flux 95 - 107 111- 98 102 - 99 A 4.85 5.50 6.11 4.14 / A d i t temperature 10.8 10.7 11.7 10.9 Rat io of qe/A 20.8 19.0 16.5 25.2 Rat io of <7h/K .087 .069 .085 .102 Rat io of w / K 309 288 282 282 Note: L\, etc. refers to norm type and first or second tr ia l . Chapter 3 - Figures 160 Figure 3.1 Thermal effects of groundwater flow in a basin with homogeneous, isotropic permeability and linear water table. Case a,b,c represent thermal disturbances for permeabilities 1 0 x 10 1 8 m 2 , 2.0 x 10 - 1 6 m 2 , 5.0 x 10" 1 6 m 2 , respectively (after Smith and Chapman, 1983) 0 0.5 1.0 Figure 3.2 Finite element mesh, various parameterization and boundary conditions for simple two dimensional synthetic examples. Case " A " , small inclusion of high log conductivity, Case d , small inclusion of high log-conductivity. Chapter 3 - Figures 162 Figure 3.3 Model misfit 5? versus Bulk Peclet number P for Case " A " , and for Case "B". Chapter 3 - Figures 163 1.0 • 0.833 1.0 0.0 0.8 0.7 0.8 0.5 0.4 0.3 0.2 0.1 0.0 0.887 -& 0.S -E 0.333 -0.187 b. 1.0 • 0.833 1.0 0.887 1 ° 0.333 -0.187 -0.9 0.8 0.7 0.8 0.4 0.2 0.0 "T 1 1 T 1 — i " 1.0 0M» 1.0 0.887 w 0.8 0.333 -0.187 0.8 0.8 0.0 1 1 1 1 0.187 0.333 0.8 0.887 0.833 1.0 mstrss Figure 3.4 ii^Tn.^ ' =1 * i o - e v * * * * * Chapter 3 - Figures 164 SYNTHETIC DATA EXAMPLE K,—7.75 K-^- - g(x,z) o 0 600 1000 1500 2000 2500 3000 DISTANCE, metre* Figure 3.5 Actual model and boundary conditions for synthetic example. Three parameter inve Aquifer location at mid-basin level. True log conductivity values of -7.75, -5.50, -20.0. Chapter 3 - Figures MK 1500-1 Figure 3.6 Hydraulic head field for aquifer at shallow depths. Chapter 3 - Figures 1600- , Figure 3.7 T e m P e m U r e fieM ! ° ' « shallow depth, . Chapter 3 - Figures 167 INVERSE SOLUTION 1600 -1 Actual: -5.75 Inverse: -5.76 1000 1500 DISTANCE, metres 2000 2500 3000 Figure 3.8 Location of aquifer as predicted by the inverse. L 2 norm functional. Aquifer at mid-basin level. Chapter 3 - Figures 168 —r 1 r 1 1 900 1000 1500 2000 2500 3000 DISTANCE, metres Figure 3.9 Locat ion of aquifer as predicted by the inverse. L2 norm functional. Aquifer at depth. Chapter 3 - Figure* 160 Figure 3.10 Groundwater flow domain with fictitious exploration program. Dots represent thermal data points, triangles are hydraulic head data points. True aquifer location shown. Chapter 3 - Figures 170 1600-1 —, , , , 1 500 1000 1500 2000 2500 3000 DISTANCE, metres Figure 3.11 Location of aquifer as predicted by the inverse. L 2 norm functional, based on hydraulic head data alone. Chapter 3 - Figures 171 I6OO-1 o H 1 1 1 1 1 \ 0 500 1000 1500 2000 2500 30*0 DISTANCE, metres Figure 3.12 Location of aquifer as predicted by the inverse. L2 norm functional, based on hydraulic heads and temperatures. Chapter 3 - Figures 172 1250 n DISTANCE, metres Figure 3.13 Log-conductivity structure from Chapter 2. _ _ _ _ Chapter 3 - Figures VtZ a. Plan LEQCMO metres Figure 3.14 Downie Slide, British Columbia, Canada, (a) Plan view of boreholes used profiled with thermal equipment, (b) Cross-section of slide indicating basal shear plane. Chapter 3 - Figures 174 Up-slope Holes 1 S8 TEMPERATURE. °C 1 X Mid-slope Holes 4 « t TEMPERATURE. °C - I — 10 —1— It I \ Bottom-slope Holes 811 814 830 832 837 840 X a. 8 Toe ***. 1-04-01C (PI) • 17 * • . - I 14 TEMPERATURE. °C - 1 — 10 TEMPERATURE. °C Figure 3.15 Thermal profiles, (a) Up-slope boreholes, (b) Mid-slope boreholes, (c) Bottom-slope bore-holes, (d) Toe boreholes. Chapter 3 - Figures 175 Figure 3.16 Downie Slide inverse. Boundary conditions and free-surface location. '500 C b * P t e r 17 u r e Geld 1SQQ DlST*NC£, ~ 30°0 2500 * B ° * » * S B * " n d i c t « ton :„ e- Con ntervay e 9 u a / s 2 o c C H A P T E R 4 Conclusions In order to solve the inverse problem in groundwater hydrology two concepts need to be considered. Firs t , a formalism or philosophy about the problem needs to be developed. In this thesis, the philosophy is adopted that a variety of models should be constructed under different norms. Only the features of the model that are consistently replicated are candidates for the 'true' model. Second, a robust inversion scheme capable of solv-ing the norms described above is needed. A groundwater inverse technique based on an opt imizat ion procedure is described. The algorithm requires no derivative computations and can be used to minimize an arbitrarily-complicated non-linear functional. The algo-r i thm produces a wide variety of acceptable models clustered about a global minimum. Each model generates data that matches observed values wi thin a noise level. Under the assumption that enough points have been sampled, the centroid of the sampled points is taken as the best estimate of the parameters. The covariances of the estimated parameters and estimated data are easily computed. Any norm can be minimized with the algorithm. Various forms of both the L\ and L2 norms are chosen depending on assumptions of the nature of noise and a priori information available. It is customary in many inversion schemes to employ kriging techniques to interpolate hydraulic heads to al l nodal points in a numerical gr id. Under noise-free conditions with adequate spacing between data points this interpolation is possible. It was observed how-ever, that on a cross-sectional problem wi th noisy data and a strong drift component, this interpolation was not possible. The difficulties are related to the inabili ty of the residual kr iging scheme to correctly identify drift in the presence of noise wi th a poorly distributed data set. 176 Conclusions 177 O n a practical problem, the technique has been applied to a sparse, inaccurate hy-draulic head data set at Downie Slide, Bri t ish Columbia . Results of a pre-inverse analysis indicate that the best that can be acheived at this site is to determine a simple model structure that explains the data. Whi le more complicted structures may exist that explain the data, the simplest structure may be the most geologically unique. Four parameters are determined; the free-surface position of the water table and three boundary conditions for the domain. Further simulations using a theoretical data set wi th assumed properties similar to that of Downie Slide show that with noise free data, and an adequate spacing between points it is possible to interpolate an unbiased estimate of hydraulic head data at all nodes in the equivalent discretized domain. When the inverse technique is applied, the domain's conductivity structure is correctly identified when enough prior log-conductivity information is available. The implications for Downie Slide are that in order to construct anything but a simple hydrogeologic model, accurate field measurements of hydraulic head are required, as well as well-defined estimates of hydraulic conductivity, a better spacing between measurements, and adequate knowledge of the boundary conditions. One way of overcoming data limitations in an inverse problem is to introduce a sec-ond data set for another potential and apply simultaneous inversion. The joint inversion method seeks to use data from a number of different measurements to improve the reso-lut ion of parameters which are in common to one or more functional relationships. One such data set is subsurface temperature, which is sensitive to variations in hydraulic con-duct ivi ty in more permeable systems. T h e nature of an advective disturbance has a significant impact on the proposed method. In low-permeability systems, where the thermal regime is conductive, inver-sion of the temperature field to construct hydraulic conductivity models wi l l not yield good results; the temperature field being far too insensitive to hydraulic conductivity vari-ations. Similarly, at high fluid velocities variations in hydraulic conductivity w i l l not effect Conclusions 178 temperature patterns. The basic idea behind the joint inversion scheme is simple; temperature measurements are taken in piezometers or boreholes along with hydraulic head measurements. Many more temperatures than hydraulic head measurements can generally be taken in most situations. This feature increases the amount of temperature data compared to hydraulic heads that can be gathered at most sites. In addition, temperature measuring equipment is simple to construct and use in the field. The simultaneous inverse scheme is tested on a variety of synthetic examples. It can be concluded from studies of simple hydraulic conductivity structures that the joint inverse scheme is sensitive to the Bulk Peclet number of the system, which is a measure of the level of advective disturbance. In other examples the ability of the simultaneous inverse scheme to reconstruct a hypothetical layered medium is examined. One of the common problems in groundwater hydrology is determining the location and hydraulic properties of aquifers in a flow system. This problem is particularly difficult if exploration boreholes do not penetrate the aquifer or where boundary conditions are uncertain. W i t h a limited data base in hydraulic head the inverse indicates a smearing of the hydraulic conductivity in zones below measurement points. In the synthetic examples chosen, the inclusion of thermal data correctly identified the recharge rate and the location and hydraulic conductivity of the aquifer. In another series of simulations the lower boundary to a flow system is determined wi th the joint inverse. A d d i t i o n a l runs on a synthetic, heterogenous medium have been carried out. W i t h a good temperature data base, thermal properties can be properly resolved. However, this stochastic problem (with certain boundary conditions) the addition of thermal data did not condit ion the inverse to a greater degree than accomplished wi th the addition of prior information on log-conductivity. The benefits of including thermal data and applying a jo int inversion can be substantial when considering the more realistic problem of uncer-ta in boundary conditions. A simulation is carried out wi th unspecified hydraulic head Conclusions 179 boundaries, an uncertain recharge rate, uncertain basal heat flux and an uncertain ther-mal conductivi ty. Good agreement between actual and computed values is indicated. This same problem based on hydraulic head data alone or with a large base in prior hydraulic conductivi ty data would be difficult to solve uniquely. The simultaneous inverse is also applied to a data set collected at the Downie Slide. Results on an inverse based on hydraulic head data alone show a simple model of four parameters: the free-surface position of the water table, and three boundary conditions for the domain. Unfortunately, with a homogeneous hydraulic conductivity, all that can be determined from the inverse are ratios of flux to hydraulic conductivity. B y including thermal data, the value of hydraulic conductivity can be determined at this site. Some of the model parameters (basal heat flux, thermal conductivity, specified head boundaries) are not resolved well by the joint scheme. None-the-less the constructed models do offer valuable insight into the hydrogeology of the field site. The constructed models persistently show a hydraulic conductivity value of about 1 x 1 0 - 7 m/sec, which is consistent wi th previous estimates of hydraulic conductivity at the site. A basal layer log conductivity is determined to be -8.6, which indicates that a low hydraulic conductivity base in the flow system is inferred. A further comparison with the hydraulic head inverse results show good agreement between the two inverses for the hydraulic properties. Recharge rates are estimated to be about 7 x 10~ 9 m/sec; very close to the estimated value of 1 x 10~ 8 m/sec based on meteorlogical studies and baseflow calculations. Basal heat flux values are found to be in the range of 114 to 100 m W / m 2 , which is reasonable based on heat flow studies south of the project area. Thermal conductivity values are determined to be about 5.0 W / m ° C sec; a value that indicates a high silica content in the rocks, more typical of quartzit ic gneisses than schists.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Simultaneous inversion of thermal and hydrogeologic...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Simultaneous inversion of thermal and hydrogeologic data Woodbury, Allan David 1987
pdf
Page Metadata
Item Metadata
Title | Simultaneous inversion of thermal and hydrogeologic data |
Creator |
Woodbury, Allan David |
Publisher | University of British Columbia |
Date Issued | 1987 |
Description | The question that is addressed in this thesis is: can a simultaneous inverse scheme involving thermal and hydrologic data resolve hydrologic model parameters to a better degree than hydrologic data alone? The first chapter sets the framework for this question by first reviewing linear and non-linear inverse problems and then illustrating the advantages of a simultaneous inverse of two different data sets through the use of a simple example. It is the goal of Chapter 2 to examine current methodologies for stating and solving the inverse problem. A review of the maximum likelihood approach is presented, and a construction formalism is adopted by introducing a series of objective functionals (norms) which are minimized to yield a variety of possible models. The inverse is carried out using a modification of a constrained simplex procedure. The algorithm requires no derivative computations and can be used to minimize an arbitrarily complicated non-linear functional, subject to non-linear inequality constraints. The algorithm produces a wide variety of acceptable models clustered about a global minimum, each of which generates data that match observed values. The inverse technique is demonstrated on a series of one and two-dimensional synthetic data sets, and on a hydraulic head data set from Downie Slide, British Columbia, Canada. At this site, four parameters are determined; the free-surface position of the water table and three boundary conditions for the domain. Further simulations using a theoretical data set with assumed properties similar to that of Downie Slide show that with noise free data, and an adequate spacing between points it is possible to interpolate an unbiased estimate of hydraulic head data at all nodes in the equivalent discretized domain. When the inverse technique is applied, the domain's conductivity structure is correctly identified when enough prior log-conductivity information is available. The implications for Downie Slide are that in order to construct anything but a simple hydrogeologic model, accurate field measurements of hydraulic head are required, as well as well-defined estimates of hydraulic conductivity, a better spacing between measurements, and adequate knowledge of the boundary conditions. Chapter 3 is devoted to developing the idea of a joint inversion scheme involving both thermal and hydrologic data. One way of overcoming data limitations (sparse hydraulic head or few hydraulic conductivity estimates) in an inverse problem is to introduce an independently collected data set and apply simultaneous or joint inversion. The joint inversion method uses data from a number of different measurements to improve the resolution of parameters which are in common to one or more functional relationships. One such data set is subsurface temperature, which is sensitive to variations in hydraulic conductivity. In Chapter 3, the basic concepts of heat and fluid transfer in porous media with emphasis on forced convective effects are reviewed, followed by inversion of theoretical data and a re-investigation of the hydrogeology of Downie Slide, augmented with thermal data and a simultaneous inverse. Additional runs on a heterogenous medium presented in Chapter 2 are carried out. With a good temperature data base, thermal properties can be properly resolved. However, in this stochastic problem the addition of thermal data did not condition .the inverse to a greater degree than accomplished with the addition of prior information on log-conductivity. The benefits of including thermal data and applying a joint inversion can be substantial when considering the more realistic problem of uncertain boundary conditions. The simultaneous inverse is also applied to the Downie Slide data set examined in Chapter 2. Unfortunately, with a homogeneous hydraulic conductivity, all that can be determined from a hydraulic head inverse are ratios of flux to hydraulic conductivity. By including thermal data, the value of hydraulic conductivity can be determined at this site. Some of the model parameters (basal heat flux, thermal conductivity, specified head boundaries) are not resolved well by the joint scheme. None theless the constructed models do offer valuable insight into the hydrogeology of the field site. The constructed models persistently show a hydraulic conductivity value of about 1 x 10⁻⁷ m/sec, which is consistent with previous estimates of hydraulic conductivity at the site. A further comparison with the inverse results in Chapter 2 show good agreement between the two inverses for the hydraulic properties. |
Subject |
Hydrogeology -- Data processing |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-08-23 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052355 |
URI | http://hdl.handle.net/2429/27652 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geological Sciences |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1987_A1 W66_2.pdf [ 7.98MB ]
- Metadata
- JSON: 831-1.0052355.json
- JSON-LD: 831-1.0052355-ld.json
- RDF/XML (Pretty): 831-1.0052355-rdf.xml
- RDF/JSON: 831-1.0052355-rdf.json
- Turtle: 831-1.0052355-turtle.txt
- N-Triples: 831-1.0052355-rdf-ntriples.txt
- Original Record: 831-1.0052355-source.json
- Full Text
- 831-1.0052355-fulltext.txt
- Citation
- 831-1.0052355.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0052355/manifest