EFFICIENT AND RESPONSIBLE USE OF PRIOR INFORMATION ANDJOINT PARAMETER ESTIMATIONFOR ESTIMATING PARAMETERS OF GROUNDWATER FLOW MODELSbyRICHARD A. WEISSB.Sc., The Colorado School of Mines, 1986M.Sc., The Pennsylvania State University, 1989A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE FACULTY OF GRADUATE STUDIESDepartment of Geological SciencesWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIAApril 1994© Richard Alfred Weiss, 1994In presenting this thesis in partial fulfUrnent of the requirements for an advanceddegree at the University of British Columbia, I agree that the Ubrary shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)Department of C€.0(LtqJ SctLRJosThe University of British ColumbiaVancouver, CanadaDate_________________DE.6 (2188)AbstractParameter estimation problems for groundwater flow systems are often non-identifiable or unstable using only hydraulic head data. Prior information onparameters or joint parameter estimation including tracer concentrations may be usedto stabilize the parameter set. Response surfaces and multiparameter confidenceregions are used to identify the most efficient and responsible use of prior informationand joint data sets for the purpose of parameter estimation.Prior information on some of the parameters are used to stabilize the parameterset for parameter estimation. Efficient use of prior information involves identifyingthose parameters for which prior information will stabilize the model parameter set themost. Responsible use of prior information involves identifying how errors in the priorinformation will influence the parameter estimates. The most responsible parametersfor prior information are those parameters for which errors in the prior informationhave the least influence on the final parameter estimates. Guidelines are developedfor efficient and responsible use of prior information in parameter estimation based onan analysis of the parameter space using response surfaces and eigenspacedecomposition.Joint parameter estimation is used when more than one data set is availableto estimate the model parameters. Response surfaces and confidence regions areused to show how multiple data sets reduce parameter uncertainty. The value offuture data in reducing the uncertainty of parameter estimates is explored. Forweighting data sets in joint parameter estimation, three criteria based on parameterspace analysis are proposed. These three criteria are evaluated and compared to themore traditional weighting method based on an analysis of data residuals.A groundwater model for the San Juan Basin, New Mexico, is constructed andcalibrated using the methods developed in this thesis. Hydraulic head data, 14Cconcentration data, and prior information on the model parameters is used to calibratethe model in an efficient and responsible manner. The model is calibrated in twostages. In the first stage, prior information on the hydraulic conductivity parametersfor the lower model layers was found to be both efficient and responsible in stabilizingthe parameter set. In the second stage, the parameter estimates and uncertaintiesbased on the four weighting criteria were compared.IIITable of ContentsAbstractTable of ContentsList of TablesList of Figures.AcknowledgementsIIivixxiixv1. INTRODUCTION1.1 Need for parameter estimation . .1.2 lIl-posedness of the inverse problem.1 .3 Uniqueness, identifiability, and stability1 .4 Methods for controlling instability .1.4.1 Parameter dimension1.4.2 Prior information1.4.3 Additional data1 .5 Outline of thesisissues136810101214172. REVIEW OF PARAMETER ESTIMATION2.1 Single state parameter estimation2.1.1 Parameterization2.1.2 Model construction2.1.3 Direct and indirect techniques2.1.4 Formulation of inverse problem2.1.4.1 Non-linear regression2.1 .4.2 Bayesian estimation2.1.4.3 Maximum likelihood estimation2.1.4.4 Stochastic formulation2.2 Joint parameter estimation2.3 Parameter space analysis 341919202224252627292930iv3. PHILOSOPHICAL AND NUMERICAL ISSUES IN THE INVERSE PROBLEM 373.1 Philosophical issues 373.1.1 Parameterization 383.1.2 Method of parameter estimation 393.2 Computational aspects of non-linear regression 393.2.1 Solution to the forward problem 403.2.2 Solution to the inverse problem 433.2.2.1 Calculation of the sensitivity matrix 463.2.2.2 Modifications to the basic Gauss-Newton method . . . 514. CONCEPTS OF PARAMETER SPACE ANALYSIS 544.1 Parameter space 544.2 Response surfaces 554.3 Relationship between response surfaces and confidence regions . . . 584.4 Linearized confidence regions 604.5 Examples of response surfaces 634.6 Relative contribution to the coefficient of variation 654.7 Summary 665. PARAMETER SPACE ANALYSIS FOR EFFICIENT AND RESPONSIBLE USEOF PRIOR IN FORMATION 715.1 Synthetic model 745.1.1 Parameter space T1-R 755.3 Efficient use of prior information . 765.3 Effect of errors in prior information . 795.3.1 Consequences of errors in prior estimates 805.3.2 Conceptual relationship between confidence region and errorratio 845.3.3 Calculation of error ratios from confidence region . 87V5.3.4 Using error ratios to identify responsible parameters for priorinformation 895.3.5 Influence of non-linear confidence regions 925.4 Error ratios for prior information with uncertainty 935.4.1 Example of the influence of the topology of the data responsesurface 965.5 Summary of guidelines for use of prior information 975.6 Application to a multi-parameter system 985.6.1 Parameter space of multi-parameter problem 995.6.2 Anticipated influence of prior information based on parameterspace 1055.6.3 Comparison of parameter space and modelled results 1065.6.4 Discussion of results for multi-parameter system 1125.7 Summary and conclusions 1146. PARAMETER SPACE METHODS IN JOINT PARAMETER ESTIMATION 1316.1 Parameter space analysis of multiple data sets 1326.1.1 Demonstration using a two parameter set 1336.1.2 Multiple parameter dimensions 1386.1.3 Multiple sampling periods 1406.2 Predicting the usefulness of a second data set 1416.2.1 Demonstration using the T1-T2 parameter set 1436.2.2 Discussion of worth of future data 1456.3 Weighting data sets in joint parameter estimation 1476.3.1 Analysis of data residuals 1496.3.2 Methods based on analysis of parameter space 1516.3.2.1 Maximum parameter stability 1516.3.2.2. Minimum total parameter uncertainty 153vi6.3.2.3 Minimizing longest axes of parameter confidenceregion 1546.3.2.4 Procedure for weighting by the parameter spacemethods 1556.3.3 Discussion of residual vs parameter space weighting criteria 1576.3.4 Comparison of weighting criteria 1596.3.4.1 Adaptability to changes in the model 1626.3.5 Discussion of results from four weighting criteria 1656.4 Analysis of multi-parameter system using joint data sets 1676.4.1 Analysis of parameter space based on concentration dataset 1686.4.2 Joint parameter estimation using head and concentration datasets 1736.4.3 Analysis of joint confidence region 1766.5 Summary 1817. JOINT PARAMETER ESTIMATION FOR THE SAN JUAN BASIN USINGPARAMETER SPACE METHODS 1937.1 Study area 1937.2 Geology 1947.3 Hydrogeology 1957.4 Model construction 2037.4.1 Model boundary conditions . . 2057.4.2 Hydrogeological units in the model 2077.4.3 Concentration boundary conditions and parameters 2087.4.4 Model grid .. 2117.5 Model calibration strategy 2117.6 Model calibration for 10 flow parameters 2137.6.1 Parameter estimates based on the hydraulic head data set . 213vii7.6.2 Parameter space based on head data set 2167.6.3 Parameter estimates based on 14C data set 2227.6.4 Parameter space based on 14C data set 2257.6.5 Joint parameter estimation 2277.6.6 Parameter space based on joint data set 2317.6.7 Discussion of results for 10 parameter system 2357.7 Model calibration for 6 flow parameters 2377.7.1 Parameter estimates and confidence region based on hydraulichead data set 2377.7.2 Parameter estimates and confidence region based on 14C dataset 2407.7.3 Joint parameter estimation for 6 parameter system 2437.7.4 Parameter space based on joint data set 2457.7.5 Effect of errors in other parameters 2477.8 Summary of model calibration for the San Juan Basin 2498. SUMMARY 262REFERENCES 264VIIIList of TablesTable 4.1 Axes of confidence ellipsoid 61Table 4.2 Relative contribution to total CV from each axis of the confidenceregion 66Table 5.1 Axes of confidence ellipse for T1-R parameter set 76Table 5.2 Relative contributions to the parameter CV from each axis of theconfidence region for parameter set T1-R 76Table 5.3 Consequences of error in specified head boundary value 81Table 5.4 Consequences of error in specified value of recharge 82Table 5.5 Axes of confidence region forT1-T2H parameter set 85Table 5.6 Axes of confidence region forT1-T2R parameter set 87Table 5.7 Error ratios for the T1-T2H parameter set 91Table 5.8 Error ratios for theT1-T2R parameter set 91Table 5.9 Parameter estimates and uncertainties for multiparameter problem 100Table 5.10 Axes of parameter confidence region for multiparameter problem 101Table 5.11 Relative contribution to the CV from each axis of the confidenceellipsoid 102Table 5.12 Error ratio matrix for multiparameter system 103Table 5.13 Calculated results from specifying some representative parametersat their true values 1 07Table 5.14 Calculated results from prior information on various parameters . 109Table 5.15 Calculated error ratios due to 25% error in the specified parameter 111Table 6.1 Parameter estimates and uncertainties for parameter sets A and B 135Table 6.2 Parameter estimates and uncertainties for joint parameter set usingfour different weighting criteria 160Table 6.3 Parameter estimates and uncertainties for joint data set using fourweighting criteria for the alternate model 164Table 6.4 Axes of parameter confidence region for multiparameter problembased on concentration data 170ixTable 7.12 Parameter estimates and uncertainties for joint data set using fourweighting criteria 228Table 7.13 Axes of confidence region for parameter estimates based on jointdata set at MINLEN weights 232Table 7.14 Relative contribution to the CV from each axis of the confidenceellipsoid based on the joint data set at MINLEN weights 233Table 7.15 Error ratio matrix for joint parameter estimates based on MINLENcriterion 234Table 7.16 Parameter estimates and uncertainties for joint data set based onfour weighting criteria; 6 parameter system 239Table 7.17 Axes of confidence region for parameter estimates based on headdata; 6 parameter system 240Table 7.18 Relative contribution to the CV from each axis of the confidenceregion based on head data; 6 parameter system 240Table 7.19 Axes of confidence region based on ‘4C data; 6 parameter system 242Table 7.20 Relative contribution to the CV from each axis of the confidenceellipsoid based on the 14C data; 6 parameter system 242Table 7.21 Axes of confidence region based on joint data set at MINLENweights; 6 parameter system 246Table 7.22 Error ratio matrix based on joint data set at MINLEN weights; 6parameter system 246Table 7.23 Error ratio matrix for transport parameters using parameter spacebased on joint data set at MINLEN weights; 6 parameter set . . . 248xiList of FiguresFigure 4.1 Response surface, univariate confidence interval and joint confidenceregion for parameters b1 and b2 67Figure 4.2 Construction of the axes of the confidence ellipse from the axislengths and unit vectors 68Figure 4.3 Example of a response surface for a poorly-conditioned parameterset 69Figure 4.4 Axes of confidence region showing partial coefficient of variation fromeach axis of the confidence region 70Figure 5.1 Five parameter synthetic flow system 117Figure 5.2 Response surface for T1-R parameter set using hydraulic headdata 118Figure 5.3 Response surface for T1-R parameter set using prior information onR 119Figure 5.4 Response surface for T1-R parameter set using head and priorinformation on R 120Figure 5.5 Response surface for T1-T2 parameter set using hydraulic headdata 121Figure 5.6 Effect of errors in head and recharge parameter values on estimatesofT1andT2 122Figure 5.7 Orientation of longest axis ofT1-T2H confidence ellipsoid 123Figure 5.8 Orientation of longest axis of T1-T2R confidence ellispoid 124Figure 5.9 Generic confidence ellipse to show calculation of error ratio forparameter b2 due to an error in parameter b 125Figure 5.10 Influence of shape of confidence ellipse on error ratio for (a) wellconditioned and (b) poorly-conditioned confidence regions 126Figure 5.11 Response surface for T1-T2 parameter set using incorrect priorinformation on T2 127xiiTable 6.5 Relative contribution to the CV from each axis of the confidenceellipsoid using the concentration data 171Table 6.5 Error ratio matrix for multiparameter problem based on concentrationdata 172Table 6.7 Parameter uncertainties for joint data set using four weighting criteriafor multiparameter problem 175Table 6.8 Error ratio matrix for transport parameters using joint data set . . 178Table 6.9 Error ratio matrix for the multiparameter problem using the joint dataset 180Table 7.1 Hydraulic head data used for calibrating the San Juan Basinmodel 197Table 7.2 The ‘4C data used to calibrate the San Juan Basin model 200Table 7.3 Prior information on the parameters in the San Juan Basin model 204Table 7.4 Parameter estimates and uncertainties based on hydraulic headdata 214Table 7.5 Fit of observed and simulated head data at the parameter estimatesbased on the head data 217Table 7.6 Axes of confidence region for parameter estimates based on headdata set 219Table 7.7 Relative contribution to the CV from each axis of the confidenceellipsoid based on the head data set 221Table 7.8 Parameter estimates and uncertainties based on 14C data set . . 223Table 7.9 Fit of observed and simulated 14C data using parameter estimatesbased on ‘4C data set 224Table 7.10 Axes of confidence region for parameter estimates based ondata set 225Table 7.11 Relative contribution to the CV from each axis of the confidenceregion based on 14C data set 226xFigure 5.12 Response surface for T1-R parameter set using head and incorrectprior information on T1 128Figure 5.13 Response surface for T1-T2 parameter set using head and incorrectprior information on parameter T2 129Figure 5.14 Multi-parameter synthetic flow system showing arrangement oftransmissivity and recharge zones 130Figure 6.1 Tracer concentration distribution for 5 parameter flow system at 500days 184Figure 6.2 Response surfaces for parameter set A: (a) head data; (b)concentration data; (c) joint data set 185Figure 6.3 Response surfaces for parameter set B: (a) head data; (b)concentration data; (c) joint data set 186Figure 6.4 Response surfaces for T1-T2 parameter set: (a) synthetic futureconcentration data; (b) head and synthetic future concentration data;(c) actual concentration data; (d) head and actual concentrationdata 187Figure 6.5 Two confidence ellipses with the same volume and orientation, butdifferent stabilities and parameter uncertainties 188Figure 6.6 Response surfaces forT1-T2 parameter set using joint data sets: (a)RESID weights; (b) MINCN weights; (c) MINVOL weights; (d)MINLEN weights 189Figure 6.7 Response surfaces for alternate model: (a) head data; (b)concentration data 190Figure 6.8 Response surfaces for alternate model using joint data sets: (a)RESID weights; (b) MINCN weights; (c) MINVOL weights 191Figure 6.9 Concentration distribution for multi-parameter flow system at 5000days 192Figure 7.1 Plan view of the San Juan Basin showing outcrops of major unitsand the location of the modelled cross section 251XIIIFigure 7.2 Generalized north-south cross section of the San Juan Basin alongmodelled cross section 252Figure 7.3 Location of wells with hydraulic head data in R9W, Ri OW and Ri 1Wfrom Chaco Canyon to the San Juan river 253Figure 7.4 Location of midpoint of screened interval for wells in head data setalong modelled cross section 254Figure 7.5 Location of wells with ‘4C data in R9W, R1OW and R11W fromChaco Canyon to the San Juan river 255Figure 7.6 Location of midpoint of screened interval for wells with 14C data alongmodelled cross section 256Figure 7.7 Conceptual cross section showing layering and model boundaryconditions 257Figure 7.8 Head distribution based on final parameter estimates using thehydraulic head data data set 258Figure 7.9 The 14C distribution based on final parameter estimates from 14C dataset 259Figure 7.10 Head distribution based on final parameter estimates from joint dataset 260Figure 7.11 The 14C distribution based on final parameter estimates from jointdata set 261xivAcknowledgementsFrom the moment I arrived at the University of British Columbia, Dr. LeslieSmith guided the way toward the completion of this thesis. He has been anoutstanding research supervisor, devoting an extraordinary amount of attention towardmy work. I would like to thank him for his time, energy, and most of all hisunderstanding. This research has been supported by University Graduate Fellowshipsfrom the University of British Columbia, as well as grants from the Natural Sciencesand Engineering Research Council of Canada (NSERC). I would like to thank MartinShute for the use of his unpublished 14C data for the San Juan Basin. I would alsolike to thank Tom Clemo for his willingness to solve any computer problems thatarose. Finally, I would like to thank my wife, Rosanna Weiss, who supported,encouraged and helped me continuously throughout this project.xvCHAPTER 1. INTRODUCTIONMathematical models are commonly used to simulate fluid flow and solutetransport in subsurface environments. Site specific models are constructed in orderto better understand the subsurface environment or to predict futurevaIues ofhydraulic head and mass concentration. The simulated model results will be accurateonly if the model is a reasonable representation of the true flow and transport system.While the mathematical and computational aspects of groundwater models are welldeveloped, the question of how to choose appropriate parameter values for a specificflow and transport system is problematic. The process of choosing appropriateparameter values for a model is termed parameter estimation.Groundwater flow can be described mathematically by a partial differentialequation and associated boundary conditions. Models for groundwater flow requireparameter values at every point in the model domain to solve for the dependentvariable, the hydraulic head distribution in space and time. These parameters includehydraulic conductivity, storage coefficient, boundary conditions and initial conditions.It is impossible to directly identify a parameter such as hydraulic conductivity at everypoint in the domain being modelled, so some method of indirectly identifying theparameters must be found. The inverse method of parameter estimation usesobservations of the dependent variable (hydraulic head) to estimate the parametersof the model.1Hydraulic head data have traditionally been used to estimate model parameters.However, the inverse problem using only head data is often ill-posed [Emsellen andde Marshy, 1971]. The ill-posed nature of the inverse problem leads to unstable anduncertain parameter estimates. Independent information (prior information) onparameter values, such as hydraulic conductivity estimates from pump tests, havebeen used to reduce the problems associated withthe ill-posed nature of the inverseproblem [Cooley, 1982]. Another way of reducing the problems due to the ill-posednature of the inverse problem is to use an additional data set through joint parameterestimation. Environmental tracers (isotopes) are present throughout manygroundwater systems, and the tracer data may be used in conjunction with hydraulichead data to provide improved parameter estimates. For large scale flow systems,environmental tracer data is a practical type of mass concentration data for use inestimating flow system parameters.The work in this thesis was motivated by an attempt to understand how priorinformation and joint parameter estimation stabilize parameter estimates and reduceparameter uncertainty. In numerical experiments using a variety of models, it wasnoted that in some cases prior information significantly stabilized the parameterestimates, while in other cases prior information did not stabilize the parameterestimates at all. Similarly, for some models joint estimation stabilized parameterestimates significantly over single state parameter estimation, while for other modelsjoint estimation did not stabilize parameter estimates at all. The work in this thesis has2focused on understanding why the above effects occur, and developing methods todetermine how prior information and joint parameter estimation can be used in anefficient and responsible manner.The most efficient parameters for prior information are those parameters forwhich prior information wiIlstabilize the parameter estimates the most. Responsibleuse of prior information involves identifying how errors in the prior information willinfluence the parameter estimates. The most responsible parameters for priorinformation are those parameters for which errors in the prior information have theleast influence on the final parameter estimates. For joint parameter estimation,parameter space analysis is used to determine the information available about themodel parameters from each of the data sets. Several criteria are developed in thisthesis to determine how to combine the information from each data set to produce thebest parameter estimates.1.1 Need for parameter estimationSteady state groundwater flow is governed by the equation:V-(K.Vh)--q=O onR (1.1)subject to the boundary condition:3-K• Vh n = Q on (1.2a)h=H onF2 (1.2b)where R is the spatial domain of the flow system, F is the boundary of the flowsystem, V is the gradient operator, K is the hydraulic conductivity tensor, h is thehydraulic head, q represents internal sources, n is the unit vector normal to theboundary, Q is the prescribed boundary flux on F1, and H is the prescribed boundaryhead on F2.The solution of the flow equation requires the knowledge of the parameters K,q, H, and 0 over the entire flow domain and on its boundary. The output from themodel are the simulated values of h throughout the domain of the model. In practice,the flow equation is often solved by numerical methods, and the parameters take ondiscrete values. These discrete parameter values will be referred to as the modelparameter values.The most straightforward way of constructing a model of a flow system is toobtain measurements of all the system parameters, use the measured parametervalues as the model parameter values, and directly calculate the output head values.There are several problems with this approach. First, due to the heterogeneity of themedia, measurements of the parameters would be required at every point in the flowdomain. It is prohibitively difficult and expensive to measure all parameters at every4point in the flow domain. Parameters such as hydraulic conductivity are expensive tomeasure in the field, and in practice measurements are relatively scarce. Second, thescale of the measured parameters is often quite different from the scale of the model.A pump test measures hydraulic conductivity at a scale that is often much smallerthan the scale of the model. These problems necessitate an alternate method ofconstructing a flow and transport model for a specific field site.A more practical approach to constructing a model is the process of modelcalibration. One way of determining if a model is a good representation of the flowsystem is to compare the output head values from the numerical model to a set offield data. If the simulated head values match the field data closely, then the modelis deemed to be a good representation of the system (at least for the purposes ofsimulating the model output). If the simulated head values are a poor match to thefield data, then the model may be a poor representation of the system. Modelcalibration is the process of constructing a model so that the simulated data closelymatch the field data.The traditional method of model calibration is manual trial and error. The valuesof the model parameters are adjusted manually until the simulated data match theobserved data. When an adequate fit between the simulated and observed data isobtained, the model is deemed calibrated. The parameter values of the calibratedmodel are the model parameter estimates. The trial and error approach is often very5helpful in understanding how a flow system behaves. However, the trial and errorprocess is often time consuming, frustrating, and subjective. It is difficult to determinewhether the parameter estimates are the best estimates, or whether other, equallygood or better sets of parameter estimates exist. The quality of the model and themodel parameter estimates are difficult to evaluate using the trial and error approach.Formalized parameter estimation methods offer a way to automate the trial anderror model calibration process. The formalized parameter estimation schemes presenta framework in which the most likely parameter estimates can be determined, and thereliability of these parameter estimates can be quantified. The calibration process isautomated, yet the modeller still retains control over the calibration process.Theparameter estimation methods are based on the solution of the inverse problem, whichis the solution of equations (1.1) and (1.2) for the parameters (independent variables),given the values of the output variable h (the dependent variable). Chapter 2 providesan overview of many of the parameter estimation methods.1.2 IlI-posedness of the inverse problem.The inverse problem is the basis for formalized parameter estimation, but it isoften ill-posed. Due to this ill-posed nature, strict conditions must be met in order toobtain a meaningful solution. For a groundwater flow problem, these conditions arethat error free values of the hydraulic head must be known everywhere in the flow6domain, and either the hydraulic conductivity or the flux must be prescribed at onepoint along each streamline [Emsellen and de Marshy, 1971]. In addition, the modelmust be a true representation of the real system.In practice, the data are generally sparse, available only at various pointlocations throughout the model domain. For a variety of reasons, the data are also inerror with respect to the model [Cooley, 1979]. Examples of data errors are: (1)measurement errors due to instrument error; (2) errors in head data due to errors insurveyed elevation; (3) areal models assume that the data is averaged over thevertical, but the measurements may not be taken over the entire vertical interval; (4)in cross section models, there is some uncertainty in the depth of measurement whenthe screened interval of the well is long; (5) parameter variations smaller than themodeled scale may cause fluctuations in data values that are not accounted for in themodel; (6) data may be influenced by a transient flow system, yet the model may beat steady state. All of these factors contribute to errors in the data with respect to themodel, and the magnitude of these errors is often unknown.In the field, parameters such as hydraulic conductivity vary continuously inspace throughout the flow system. Numerical models of the flow system, either finitedifference or finite element models, divide the flow system into elements. Thecontinuous variation of the field parameter is replaced by a set of discrete modelparameters. The model cannot represent the true variation in the field parameter, and7is therefore an approximation to the true flow system. For groundwater problems, theabove factors lead to an ill-posed inverse problem. The ill-posed nature of the inverseproblem often makes the inverse problem difficult to solve.1.3 Uniqueness, identifiability, and stability issuesThree major issues result from the ill-posed nature of the inverse problem:uniqueness, identifiability, and stability. A detailed mathematical analysis of theseissues is contained in Carrera and Neuman [1986b]. We present these issues in aconceptual context.A solution to the inverse problem is said to be unique if the set of estimatedparameters is the only set which satisfies the conditions to be a solution [Carrera,1988]. Two types of non-uniqueness have been recognized [Tarantola and Valette,1982]. One type is due simply to the paucity of data, that is, more parameters are tobe identified than data available. The second type of non-uniqueness is due to thestructure of the physical model. Multiple minima in the objective function surface area symptom of non-uniqueness due to model structure. Non-uniqueness due to modelstructure may be relatively rare [Carrera, 1988].A model parameter is non-identifiable if the model output is not sensitive to theparameter [Chavent, 1979]. For instance, in a one-dimension flow system with8prescribed head boundaries and a uniform hydraulic conductivity, the hydraulicconductivity is non-identifiable using head data. The value of hydraulic conductivityhas no influence on the head distribution. If a parameter is non-identifiable, then thesolution to the inverse problem is non-unique. However, a non-unique solution mayoccur even if all parameters are identifiable.A group of parameters may also be non-identifiable. For example, when allmodel parameters are estimated using only head data, the set of model parametersare non-identifiable [Cooley and Naff, 1990]. The model output is not sensitive to theentire set of model parameters, considered jointly. A suite of different combinationsof parameter values result in identical model outputs. Because of the non-identifiabilityof the entire set of model parameters, all model parameters are very seldomestimated. It is common practice to estimate only a subset of the model parameters,usually the internal model parameters. The model boundary conditions are commonlyspecified.The major difficulty associated with parameter estimation arises from theinstability of the parameter set being estimated. When a parameter set is unstable,small errors in the data can lead to large errors in the estimated parameters. Sinceerrors in the data are unavoidable, large errors in estimated parameters are commonin unstable problems. Unstable parameter sets also lead to large uncertainties in theestimated parameters. Unstable parameter sets are characterized by a large, nearly9flat region in the response surface (see Chapter 4). A side effect of parameterinstability is that it slows the convergence rate of the iterative algorithms used forparameter estimation. Instability often leads to iterative algorithms that stop at pointsthat are not the true minimum, leading to the ‘(false) assumption that multiple localminima exist. Thus, instability is often confused with non-uniqueness. Someresearchers have concluded that instability in the inverse problem is such a majorproblem that sensible parameter estimates are unattainable in many instances[Yakowitz and Duckstein, 1980].1.4 Methods for controlling instabilitySince parameter instability is a major issue in parameter estimation forgroundwater flow problems, many researchers have focused on methods of controllingthe parameter instability. Three approaches are possible: (1) reduction in parameterdimension; (2) incorporation of prior information on parameter values; and (3)collection of additional data.1.4.1 Parameter dimensionThe parameter dimension is the number of parameters to be estimated. In anumerical model, an independent model parametermay be assigned to each elementfor each type of parameter [eg. Shah et al., 1978]. In this case, there will generally be10more parameters than data, and the resulting parameter estimates will be non-unique.To reduce the parameter dimension, elements may be grouped in zones, and a modelparameter assigned to each zone [eg. Cooley, 1977]. Zones for one type ofparameter, such as hydraulic conductivity, may not correspond to zones for anothertype, such as recharge.As the number of zones increases, the fit of the simulated data to the observeddata becomes better, while the uncertainty of the parameter estimates increases. Asthe number of zones decreases, the fit of the simulated data to the observed databecomes worse, but the parameter estimates have lower uncertainties. Severalmethods have been developed to determine the optimum number and arrangementof the zones (Shah et a!., 1978; Yeh and Yoon, 1981; Sun and Yeh, 1985: Carreraand Neuman, 1986a; Cooley et al., 1986). All methods attempt to minimize theuncertainty of the parameter estimates while maximizing the fit of the data to themodel. The methods select the simplest model compatible with the available data.Reducing the number of zones (parameter dimension) may not always help toreduce the problems associated with the inverse problem. For the case of estimatingall model parameters, the inverse problem is non-identifiable regardless of the numberof parameter zones. In other cases, the geology of the site may impose a givennumber and arrangement of zones. It may not be possible to reduce the number ofzones and still retain a credible model of the site.111.4.2 Prior informationIndependent information on the estimated parameters (called prior information)has been recognized as valuable in stabilizing the inverse problem. The priorinformation can either have a specified reliability [Neuman and Yakowitz, 1979] or anunknown reliability [Cooley, 1982]. For the case of estimating all the parameters fora groundwater flow model, the inverse problem becomes identifiable in the presenceof prior information [Carrera and Neuman, 1986b]. Although identifiable, the inverseproblem may still be very unstable.Though prior information may be valuable in stabilizing the inverse problem, acautious use of prior information may be advisable. Prior information may notsignificantly stabilize an inverse problem. There is little advantage to be gained byobtaining prior information about a parameter which will not be effective in stabilizingthe model parameter set. In this thesis, guidelines are developed to identify thoseparameters for which prior information will most efficiently stabilize the parameter set.Another concern when using prior information is that it may not berepresentative of the model parameter values. The prior information may be in erroror biased with respect to the model for a number of reasons. First, prior informationmay be obtained at a scale which is different from the scale of the model. As anexample, suppose the prior information about hydraulic conductivity for a model zone12is obtained using a slug test. The slug test samples smaller volume of the subsurfacethan the model zone, and the slug test may miss some large scale features that themodel zone includes. The hydraulic conductivity obtained from the slug test wouldmost likely be different than the hydraulic conductivity at the scale of the model. Ifprior information about all hydraulic conductivity zones were obtained from slug tests,the prior information would most likely be consistently biased with respect to the modelhydraulic conductivity values. Second, prior information on parameters such ashydraulic conductivity is often extrapolated from measurements taken outside of themodelled region. This extrapolation introduces the possibility of errors in the priorinformation. Third, prior information may be obtained in parts of hydrogeological unitsthat are not representative of the entire unit. As an example, in low permeability units,hydraulic conductivity measurements are often obtained in the higher permeabilitylenses within the units. Also, if the prior information is obtained from tests near theoutcrop of a hydrogeological unit, it may not represent the parameter values for thesame unit at depth. All of these factors contribute to the fact that the prior informationmay be unrepresentative of the model parameter values.The influence of unrepresentative or biased prior information on the final modelparameter estimates is explored in this thesis. The biased prior information influencesnot only those parameters with prior information, but the other model parameterestimates as well. In fact, small errors in prior estimates may lead to greatly magnifiederrors in the final estimates for other parameters. In this thesis, guidelines are13developed to identify those parameters for which errors in the prior information leadto the smallest possible errors in the final parameter estimates.1.4.3 Additional dataAdditional data on the dependent variable may be used to stabilize the inverseproblem. This additional data may take the form of additional hydraulic head data, orit may be a new data type, such as mass concentration data. If the additional data aremore hydraulic head data, they may or may not assist in stabilizing the inverseproblem. When estimating a subset of the model parameters, the additional head datamay help stabilize the parameter estimates. However, when estimating all modelparameters, the additional head data will not stabilize the inverse problem. In thiscase, the inverse problem is non-identifiable regardless of the number of head dataused.Additional data may take the form of mass concentration data. The massconcentration data will yield information on the flow parameters, since theconcentration is dependent on the fluid velocities, which are a function of the flowparameters. For this reason, the mass concentration data should assist in stabilizingthe parameter estimates for the flow parameters, even when estimating all model flowparameters. The use of head and concentration data together falls under the conceptof joint parameter estimation.14Solute transport of a single decaying nonreactive species is governed by theadvection-dispersion equation:V •9 (D•Vc - Vc) - = e. on R (1.3)subject to the initial conditions:c = c on R at t—_t0 (1.4a)and the boundary conditions:(0 DVc - OVc)•n = g1 on (1.4b)c=g2 onF4 (1.4c)where D is the dispersion tensor, c is the mass concentration of the solute, V is thefluid velocity, 2. is the decay coefficient, 0 is the porosity, c0 is the initial concentration,g1 is the specified mass flux along F3, and g2 is the specified concentration along F4.The components of the fluid velocity are calculated from the solution to the flowequation byVi= EJ13[——---—] (1.5)The components of the dispersion coefficient are given byD11 T”i,j + (cL—cT).YL + D,,611 (1 .6)where cL is the longitudinal dispersivity, cia- is the transverse dispersivity, V is themagnitude of the fluid velocity,°d is the molecular diffusion coefficient, and ö, is the15Kronecker delta.The solution to the solute transport equations requires the knowledge ofadditional parameters, such as cxi, cIT, 0, and the initial and boundary conditions onmass concentration. These additional parameters may either be estimated during thesolution of the inverse problem, or independent information on these additionalparameters could be obtained.The concentration data may eliminate the non-identifiability of the parameterset and stabilize the parameter estimates. In some instances, concentration data willstabilize the parameter estimates and reduce parameter uncertainty significantly, whilein other cases the concentration data may not stabilize the parameter estimates[Weiss and Smith, 1993]. Methods of determining whether the additional data willsignificantly stabilize the parameter estimates are developed in this thesis.For joint parameter estimation, relative weights for each data set are required.Three methods of weighting are developed in this thesis, based on obtaining the bestpossible estimates given the data available. Another issue involved in incorporatingconcentration data is the question of the additional transport parameters needed tosimulate the mass concentrations. If these additional parameters are specified usingprior information, the parameter values may be in error. Methods are developed in thisthesis to evaluate the influence of errors in the transport parameters on the estimates16of the flow parameters.15 Outline of thesisThe objective of the thesis research is to develop methods for stabilizingparameter estimates for groundwater flow models in an efficient and responsiblemanner. For determining the usefulness of prior information, guidelines are developedto identify those parameters for which prior information will most efficiently stabilizethe parameter set. Since the prior information may be in error with respect to themodel, guidelines are also developed to identify those parameters for which errors inthe prior information leads to the smallest possible errors in the final parameterestimates. Using prior information on the basis of these guidelines will result in anefficient and responsible method of stabilizing the inverse problem using priorinformation.For joint parameter estimation, methods are developed to determine theusefulness of mass concentration data for the purposes of stabilizing the estimatedparameter set. Several different weighting methods for the data sets in joint estimationare developed and evaluated. Methods of determining the influence of errors in thetransport parameters on the estimates of the flow parameters are also developed. Allof the above issues are explored through an examination of the model parameterspace.17The thesis is organized as follows. Chapter 2 contains an overview of previouswork in single state and joint parameter estimation, as well as parameter spaceanalysis. Chapter 3 details the methods and philosophy of parameter estimation usedin this work. An introduction to the concepts needed for parameter space analysis isprovided in Chapter 4. In Chapter 5, the parameter space analysis is used to evaluatethe use of prior information in parameter estimation. Guidelines are developed forefficient and responsible use of prior information. In Chapter 6, joint parameterestimation is analyzed in the context of parameter space analysis. A groundwater flowmodel for the San Juan basin, New Mexico, is constructed and calibrated to illustratethe methods developed in this work, and the results are presented in Chapter 7.18CHAPTER 2. REVIEW OF PARAMETER ESTIMATIONParameter estimation for groundwater flow systems has been investigated forover three decades. The majority of the research has focused on single stateparameter estimation, generally estimating flow parameters using a set of hydraulichead measurements. Less work has been done on coupled inverse problems, suchas joint flow-mass transport or joint flow-thermal transport problems. This chapterdescribes past research on the major issues involved in single state and jointparameter estimation for groundwater models, as well as past work on analysis ofmodels using parameter space methods.2.1 Single state parameter estimationA great deal of research has been done on parameter estimation for a singlestate variable in groundwater flow systems. Complete reviews may be found in Yeh[1986] and Carrera [1988]. The following is a short overview of the major issuesinvolved in parameter estimation for a single state variable. The first issue,parameterization, concerns how the true parameter structure is represented by thediscrete model parameters. The second issue, model construction, includesdetermining the correct model structure. The other issues concern estimating modelparameters, assuming a given parameterization and model structure.192.1.1 ParameterizationThe field parameters, such as hydraulic conductivity, vary continuously inspace. The boundary parameters, such as specified head, also vary continuouslyalong the boundary. This continuous variation in the field parameter is replaced by aset of discrete model parameters. Parameterization refers to the method by which thetrue parameter variation is represented by a set of discrete model parameters. Thereare several popular parameterization methods:1. Zonation: The flow region is divided into a number of subregions, or zones, and aconstant parameter value is assigned to each region [eg. Cooley, 1979]. The zonationapproach is conceptually simple, and the zonation may be based on geological orhydrogeological information. Since the number of parameters is less than the numberof data, the reliablifties of the parameter esfimates may be quantified. The majorcriticism of the zonation approach is that the zonation may be incorrect orunrepresentative of the true system. If the zonation is incorrect, then the parameterestimates for the zones may be incorrect or meaningless. The issue of determiningthe correct zonation is discussed under model construction below.2. Interpolation: This approach is generally used with finite element methods. A valueof the parameter is assigned to each node of the element, and interpolated within theelement by a local basis function [eg. Yeh and Yoon, 1981]. The elements used for20parameter interpolation need not be the same as the elements used for the headinterpolation.3. Discretization: The discrete values of the parameters at every model node or withineach element are taken as the model parameters. This approach can be used witheither zonation or interpolation, and results in a large number of model parameters[eg. Shah et a!., 1978]. The advantages of this method are that the model structuredoes not need to be determined a priori. The disadvantages are that a large numberof model parameters need to be estimated, often greater than the number of data.The inverse problem is underdetermined, and estimates of the parameter reliability aredifficult to obtain.4. Stochastic: The model parameters are viewed as random functions, characterizedby a mean and covariance. The estimated parameters are the mean, trend and thecovariance of the field parameters [eg. Kitanidis and Vomvoris, 1983].All of the above parameterization methods are valid, and each has advantagesand disadvantages. The method chosen is often determined by the type of informationrequired about the system being modelled and the purpose for which the model isbeing developed.212.1.2 Model constructionIn order to adequately calibrate a numerical model for a flow and transportsystem, three major model construction steps need to be followed. First, the modelshould include all physical processes relevant to the flow and transport system.Second, the structure of the model should resemble that of the real system. Third, thevalues assigned to the model parameters need to be similar to the true systemparameter values. Most of the work in parameter estimation for groundwater flowmodels has concentrated on the third step.The physical processes of the system are represented in the model by theequations used for the numerical model. For steady state groundwater flow, equations(1.1) and (1.2) adequately represent the physics of the flow process. For solutetransport, the advection-dispersion equation (1.3) has traditionally been used torepresent the physics of the mass transport process. However, it is recognized thatthis representation using a constant dispersivity is not generally valid in aheterogeneous porous medium [eg. Geihar and Axness, 1983].The model structure needs to be defined correctly. The type of boundaryconditions and the location of the model boundaries need to be determined. Within themodel domain, the structure of the hydraulic conductivity, recharge, porosity, anddispersivity distributions need to be defined. When discrete values of the parameters22within each element of the model grid are sought, it is not necessary to determine themodel structure before estimating the parameters. Since each element is a separatemodel parameter, the results of the parameter estimation will determine the modelstructure. More commonly, a zonation or interpolation scheme is used to parameterizethe model. The shape of the parameter zones and the location of the boundaries ofthe zones need to be determined before the parameters are estimated. In somecases, the geology of the site will greatly assist in determining the shape and locationof the internal zones. The model structure is defined by the shape and location of allmodel parameter zones.One approach for determining the correct model structure consists of selectingthe model structure which leads to the best model fit while maintaining parameterstability [Shah et al., 1978; Yeh and Yoon, 1981]. Alternative models with differentparameter structures are compared. Sun and Yeh [1985] have proposed a method forautomatically determining the model structure without having to construct alternativemodels. However, their approach does not include prior information on the modelparameters.Linear hypothesis testing can be used to compare alternative models. Cooleyeta!. [1986] sequentially compared pairs of models for a large scale flow system usinghypothesis testing. A simple model for the flow system was initially proposed, andmore complex models were compared to this simple model. Models that significantly23improved the fit of the data were accepted. More complex models were often rejected,based on the fact that the fit to the data was not improved significantly.Criteria-based approaches have been used to compare model structures.Carrera and Neuman [1986a] present four criteria, and use these criteria todiscriminate between competing models. The criteria tend to select the simplestmodels that are compatible with the data. Additional complexity can be built into themodels as the quantity and quality of the data increase.Once the model structure has been determined, the parameters of the modelneed to be estimated. The remainder of this chapter concerns estimating theparameters of the model.2.1.3 Direct and indirect techniquesThe techniques used to solve the inverse problem can be classified as eitherdirect or indirect. The direct technique treats the model parameters as unknowns, anduses the hydraulic head distribution and the spatial derivatives of this distribution tosolve for the unknown parameters directly. The solution of the direct techniquerequires that an estimate of the hydraulic head be known everywhere in the system,and that a value of the transmissivity is known on every streamline in the system.There are methods of relaxing these requirements; for instance, using a flatness24criterion [Emsellem and de Marshy, 1971]. Most of the early work with parameterestimation used the direct formulation of the inverse problem [Stallman, 1956]. Manymathematical techniques can be applied to the direct method of parameter estimation,among them are linear programming [Klelnecke, 1971]; quadratic programming [Hefezet al., 1975]; and kriging with matrix inversion [Yeh et al., 1983].The indirect technique uses an output error criterion, usually the differencebetween the measured head values and the model generated head values. The modelparameters are iteratively updated until the output error is minimized. The indirecttechnique has been more popular because of its flexibility and stability. Priorinformation on parameters can be incorporated into the parameter estimation scheme,and several data types can be used. The indirect technique is implemented in the formof minimizing an objective function, resulting in a type of least squares minimizationproblem. Mathematical programming techniques for minimizing the objective functionare often used, such as gradient search procedures [Jacquard and Jam, 1965];quadratic programming [Yeh, 1975]; Gauss-Newton method [Cooley, 1977]; conjugategradient method [Neuman, 1980]; maximum likelihood method [Carrera and Neuman,1 986a]; and the direct search method [ Woodbury et a!., 1987].2.1.4 Formulation of inverse problemThe specific aim of indirect parameter estimation techniques is to minimize the25difference between the measured and model generated hydraulic heads. However, thereal objective of parameter estimation techniques should to be to estimate meaningfulparameters. It is necessary to recognize that the hydraulic head data generallycontains errors, and the parameter estimation needs to be performed in a statisticalframework. Three primary statistical approaches- have been introduced: non-linearregression, Bayesian, and maximum likelihood methods. Each of these approachesstarts from a different philosophical background, however, they each result in a similarset of equations for parameter estimation. A fourth approach, based on the stochasticformulation, has also been proposed.2.1.4.1 Non-linear regressionNon-linear regression generally results in a weighted least squares criterion forparameter estimation. The weighted least squares criterion attempts to match themodel calculated hydraulic heads to the observed hydraulic heads in a least squaressense. For hydraulic head data and prior information on parameters, the objectivefunction to be minimized is:S [(h_hjTWh(hh*) + (pp*)TW( p* j (2.1)where h* is the vector of observed heads, h is the vector of model calculated headsat the observation points, W is a weighting matrix, T is the transpose operator, p* isthe vector of observed prior information on parameters, p is the model calculated prior26information, W is a weighting matrix, and ? is the trade-off between head and priordata. When S is minimized, the differences between the observed and calculatedheads and the prior and calculated parameter values are minimized.The main advantage ofthis approach is that it allows the use of a large bodyof literature [eg. GraybilI, 1 97] on hypothesis testing, definition of confidence regions,and model selection. The reliability of the estimated parameters can be determinedfrom the linearized covariance matrix [Yeh and Yoon, 1981]. Non-linear confidenceintervals for the estimated parameters can also be determined [ Vecchia and Cooley,1987; Cooley, 1993].2.1.4.2 Bayesian estimationThe Bayesian approach assumes that the data on hydraulic heads and priorinformation on parameters can be represented as joint probability density functions,and the resulting parameter estimates have a probability density function. The mostgeneral approach is presented by Tarantola and Valette [1982], and the approach ofGavalas et a!. [1976] is nearly equivalent. The observed hydraulic heads can bedescribed by a probability density function r(h). The prior information on parametersis described by the a priori probability density function r(p). The calculated dataobtained from the model is described by a conditional probability density function,Q(h/p), that is, the probability that the head values are h, given that the model and its27parameter values are described by p. The posterior conditional probability densityfunction of the parameter values, s(p/h), is the solution to the parameter estimationproblem. Bayes theorem states that s(p/h) = Q(h/p)r(p)/r(h), thus the probabilitydensity function of the parameter estimates can be constructed by knowing thestatistics of the prior information and the model.If the parameters and measurement errors are normally distributed and themodel is linear in its parameters, it can be shown that the posterior conditionalprobability density function s(p/h) is normally distributed as:s(plh) exp{_.[(h_h*)TCh1(h_h*) + (p_p*)TC(p_ *)]) (2.2)where Ch1 is the hydraulic head covariance matrix, and C1, is the parametercovariance matrix. The entire probability density function can be obtained. Theestimated parameters can be described by a probability density function, but it isdifficult and expensive to determine the entire distribution. The practical approach isto determine the maximum a posteriori point of the distribution. Maximum aposterioriestimation of a normal distribution is equivalent to minimizing the quadratic objectivefunction:S = (hhjTC ’( h’) + (pp)Tcl(pp.) (2.3)This objective function is identical to the objective function found using a least squarescriterion, when the inverse covariance matrices equal the weight matrices. The least28squares criterion and the Bayesian estimation criterion with normally distributed dataerrors and parameters lead to minimizing the same objective function.2.1.4.3 Maximum likelihood estimationMaximum likelihood estimation calls for choosing parameters that maximize thelikelihood function. If the data errors are normally distributed, the likelihood functionis similar to equation (2.2) above. The chief differences between maximum likelihoodestimation and Bayesian estimation are philosophical. Bayesian formulations considerthe model parameters to be random variables, while maximum likelihood formulationsconsider the model parameters to be fixed but unknown. The lack of informationcauses the parameters to be uncertain, but the model itself is deterministic. Maximumlikelihood estimation does not require that the model be capable of reproducing thetrue system exactly, because the parameters are chosen as most likely within theframework of a specific model structure. This recognition that the model is not exactleads to criteria for comparing competing models [Carrera and Neuman, 1 986aJ.2.1.4.4 Stochastic formulationStochastic formulations of the inverse problem emphasize capturing the spatialvariability of the parameters. The parameters are considered random functions, andare estimated at every point in the aquifer. A statistical model is proposed to represent29the spatial variability of the parameter within the flow system. This statistical modelcan generally be characterized by its mean, covariance, and trend. The spatialvariability of the heads are derived from the spatial variability of the parametersthrough a linearized flow equation. The parameters describing the mean, covarianceand trend of the model parameters are estimated using a maximum likelihoodapproach, conditioned on the head and prior parameter data. Kriging is applied toprovide minimum variance, unbiased, point estimates of the parameters using allavailable information. This approach has been described by Kitanidis and Vomvoris[1983] for a one dimensional flow system and Hoeksema and Kitanidis [1985] for twodimensional flow systems, while Dagan [1985] used analytical techniques forparameter estimation in an infinite medium and with statistically uniform parameters.Dagan and Rubin [1988] extend this approach to identifying different parameter typesin an unsteady flow regime.2.2 Joint parameter estimationJoint parameter estimation methods for groundwater flow systems have mainlybeen studied in the context of joint flow-mass transport systems and joint flow-thermaltransport systems. The joint parameter estimation methods were used becausemultiple data sets were available, and the joint estimates often resulted in lowerparameter uncertainties than single state parameter estimates. Joint parameterestimation methods often allowed more parameters to be estimated than equivalent30single state parameter estimation methods.In many parameter estimation studies using concentration data, the flow fieldand associated flow parameters are considered known, and only the transportparameters are estimated using the concentration data. Umari et at [1979] usedconcentration data with a linearized version of the transport equation to estimatedispersivities. No estimates of the reliabilities of these estimates were possible in theirformulation. Wagner and Gorelick [1986] estimated all transport parameters for a onedimensional system using concentration data, and included estimates of the reliabilityof the parameter estimates. Carrera and WaIters [1985] estimated the transportparameters and boundary conditions for a one dimensional tracer test.True joint parameter estimation studies estimate both the flow and transportparameters simultaneously. Strecker and Chu [1986] combine a method ofcharacteristics solute transport simulation with quadratic programming to estimatetransmissivities and dispersivities for a hypothetical aquifer. Because the observationsof head and concentrations had no error, they were able to accurately estimate theunknown parameters. Wagner and Gorelick [1987] considered optimal groundwaterquality management under parameter uncertainty in which a coupled flow-masstransport model was identified. Mishra and Parker [1989] considered parameterestimation for coupled unsaturated flow and transport. Coupled inverse problems ingeneral were discussed by Sun and Yeh [1990a], and equations and synthetic31examples for adjoint state parameter estimation were derived for both joint flow-masstransport and flow-thermal transport. Sun and Yeh [1990b] discussed identifiabilityissues for joint flow-mass transport systems. Galley et al. [1992] used an analysis ofdata residuals to determine the relative weights for the two data sets in applyingparameter estimation for joint flow-mass transport to a field problem. Xiang et a!.[1993] developed a composite L1 estimator to identify parameters in a joint flow-masstransport model, and showed the L1 estimator to be more robust in handlingobservation data containing outliers than the more traditional L2 estimators.Temperature data has been used in a joint inversion scheme to improve theparameter estimates [Woodbury and Smith, 1988; Wang et al., 1989]. Using onlyhydraulic head data, Woodbury and Smith [1988] could estimate only the ratiobetween recharge and hydraulic conductivity. The introduction of temperature dataallowed independent estimates for recharge and hydraulic conductivity. The solutionto the joint inverse problem did result in an improvement in the parameter estimates[Wang et a!., 1989], but temperature data is not readily available at many field sites.Use of temperature data for parameter estimation in groundwater flow systems isrestricted to more permeable systems with a significant component of verticalgroundwater flow. Joint geophysical-hydrological parameter estimation has also beenproposed by Copty et a!. [1993]. They first estimated the flow parameters using onlyhead and permeability data, then included seismic velocity data through Bayesianupdating. Semiempirical relationships were used to correlate the seismic velocity and32and hydraulic properties.Less formal methods of parameter estimation using environmental tracer datahave been employed. Campana and Simpson [19831 used 14C data to determineresidence times, vertical flow velocities, and average recharge rates for an aquiferusing discrete state compartmental models. Phillips et al. [1989] used hydraulic headand ‘4C to determine the hydraulic conductivity distribution in two confined aquifers inthe San Juan Basin, New Mexico. Dispersion was not considered when determiningthe 14C travel times. They used the estimates of hydraulic conductivity determinedfrom isotope travel times and the hydraulic head distribution to construct a numericalmodel of flow in the aquifers. Krabbenhoft et al. [1 990a,1 990b] used stable isotopesto estimate the discharge rates of a kettle lake to the surrounding aquifer. They usedtwo methods to determine the recharge rates, an isotope mass balance method anda calibrated quasi three dimensional flow and transport model. The recharge rate toa shallow unconfined aquifer system has been investigated using hydraulic head dataand tritium data [Robertson and Cherry, 1989]. A groundwater flow model wasconstructed using permeameter estimates of hydraulic conductivity, and rechargerates determined from the tritium distribution. The simulated hydraulic heads using thismodel compared well with the field measured hydraulic heads. The dispersivity of themedium was then calculated by comparing field measured tritium distributions to a onedimensional solution to the advection-dispersion equation. These less formal methodsuse the idea that environmental tracers provide additional information about the flow33system, but they lack a framework in which to quantify both parameter uncertainty andthe role of the environmental tracer data in reducing parameter uncertainty.2.3 Parameter space analysisThe results of most parameter estimation methods are typically a set ofparameter estimates and associated uncertainties. Analysis of the parameter spaceallows the modeller to visualize the results of parameter estimation, and develop aclearer understanding of the results. Parameter identifiability and uniqueness issuescan be clearly visualized using parameter spaces. Most parameter space analysisconsists of visualizing response surfaces (objective function surfaces) in twodimensional parameter spaceSorooshian and Arfi [1982] used two dimensional response surfaces to studythe sensitivity of model parameter estimates in rainfall-runoff models. They believedthat individual parameter sensitivity, measuring the sensitivity of one parameter whileholding the other parameters constant, was a poor method of determining parametersensitivity. They introduced two indices to measure two parameter concurrentsensitivities, based on the shape and orientation of the response surface in parameterspace.The relationship between model structure and parameter observability and34uniqueness in rainfall-runoff models was investigated by Sorooshian and Gupta[l 983].Response surface studies were used to investigate the effects of different modelstructures and different objective functions. Sorooshian and Gupta [1985] focused onevaluating a model based on whether the parameters could be identified. Theyintroduced some multiparametermethods for evaluating the identifiability of modelparameters. Reparameterization of rainfall-runoff models was suggested in order toimprove parameter identifiability.Two dimensional parameter uncertainty regions in rainfall-runoff models wereinvestigated by Kuczera [1 983a]. Kuczera [1 983b] evaluated the effects of includingindependent information on the parameters. This study also showed the effects ofdifferent kinds of data on parameter uncertainty, as well as the compatibility ofdifferent kinds of data.In groundwater modeling studies, parameter space or response surfaceanalysis is rare. Carrera and Neuman [1 986b] showed, using response surfaces, thatlog-transformation of the hydraulic conductivities generally resulted in betterconditioned parameter estimates. They also used two dimensional response surfacesto demonstrate an example of parameter non-uniqueness. Toorman eta!. [1992] usedparameter response surfaces to evaluate the influence of two data types on parameterestimates in unsaturated one-step outflow experiments. They evaluated the influenceof outflow and matric potential data on three model parameters, and concluded that35using both data types together generally resulted in more unique solutions. Parameteruncertainty was not evaluated.This thesis extends the use of response surfaces and parameter space analysisin several ways. First, most of the work using response surfaces described above wasrestribted to two dimensional parameter space, since it is easy to visualize in twodimensions. Eigenspace decomposition is used in this thesis to visualize multidimensional response surfaces. Second, response surfaces and parameter spaceanalysis are used to identify those parameters which most efficiently stabilize theparameter estimates. Third, parameter space analysis is used to determine the effectsof errors in prior information, and to identify those parameters which lead to thesmallest errors in the parameter estimates. Finally, response surfaces are used toshow how joint data sets reduce parameter uncertainty, and parameter space analysisis used to develop several different weighting criteria in joint parameter estimation.36CHAPTER 3. PHILOSOPHICAL AND NUMERICAL ISSUES IN THE INVERSEPROBLEMThis chapter covers the philosophical and computational aspects of parameterestimation used in this work. The choice of philosophy is often controlled by theobjectives of the study. Once a philosophy has been chosen, it controls the methodsand results of parameter estimation, as well as the interpretation of the results. Thetype of parameterization and type of estimation are the philosophical issuesconsidered. Computational aspects include the numerical scheme used forthe forwardsimulations of flow and transport and the modifications to the basic parameterestimation method.3.1 Philosophical issuesThis work is geared toward the application of parameter estimation to jointgroundwater flow and mass transport models for large scale flow systems. Theunderlying philosophy throughout the work is to develop responsible parameterestimates with minimum uncertainty. To accomplish these goals, a parameterestimation method that allows the calculation of estimated parameter reliabilities isrequired. The number of parameters needs to be smaller than the number of data inorder to calculate the estimated parameter reliabilities. A zonation approach toparameterization will be used, with the zonation controlled by the geology of the site.373.1.1 ParameterizationIn numerical models, the continuous spatial variation in the field parametersmust be replaced by a discrete set of model parameters. Methods of accomplishingthis parameterization were discussed in Chapter 2. In this thesis, zonation is used toparameterize the internal model parameters and a combination of zonation andinterpolation is used to parameterize the boundary conditions. This type ofparameterization is chosen because zonation allows the number of parameters to beless than the number of data, which in turn allows the reliabilities of the parameterestimates to be calculated. If the number of parameters were greater than the numberof data, parameter reliabilities could not be calculated.The major criticism of the zonation approach is that the zonation may beincorrect or unrepresentative of the true system. If the zonation does not reflect thetrue structure of the system, the parameter estimates for the zones may be incorrector meaningless. The philosophy of zonation taken in this work is that the geologicaldata can be used to define the large scale zones. Geological or hydrostratigraphicunits are mapped and used as the basis for zonation. There is often a sharp contrastin parameter values between geological units. Within each of these units, theparameter values vary continuously in space. If a trend in the parameter values withineach unit is hypothesized, then the geological unit is subdivided into zones to reflectthe trend in the parameter value. If no trend in the parameter value is hypothesized,38then the geological unit is not subdivided. Alternative models using different parameterstructures may be compared to decide if a geological unit should be subdivided intomultiple zones. Hypothesis testing [Cooley et a!., 1986] or criteria based methods[Carrera and Neuman, 1 986c] can be used to select the proper zonation within eachunit.3.1.2 Method of parameter estimationAny of the three statistical philosophies for parameter estimation introduced inChapter 2 can be used, since they all lead to similar objective functions. In this thesis,a non-linear regression approach is adopted, with a linearized approximation for theparameter reliabilities. This approach allows the calculation of parameter uncertaintiesand parameter confidence regions, which are used extensively in this work.3.2 Computational aspects of non-linear regressionIn order to estimate parameters, two problems must be solved; the forwardproblem and the inverse problem. The forward problem consists of solving equations(1.1) and (1.3) for the dependent variables, hydraulic head and tracer concentration.The inverse problem estimates the independent variables in equations (1 .1) and (1 .3)using measurements of head and concentration. The forward problem is solved usingfinite element methods. The inverse problem is solved using a modified Gauss-Newton39scheme.Three unique aspects are introduced in this thesis for the parameter estimationprocess. First, the forward transport problem is solved using the Arnoldi algorithm[Woodbury et al., 1990], which allows rapid calculation of tracer concentrations insimulations with many time steps. The Arnoldi algorithm is extended to permit thecalculation of sensitivity coefficients through the sensitivity equation method. Second,the matrices in the Gauss-Newton method are scaled using the parameter values.This scaling is performed so that the matrices are independent of the differences inparameter values, which is helpful when evaluating parameter uncertainties. Third, theGauss-Newton equations are solved using singular value decomposition [Press et a!.,1986]. The singular values and vectors are used to allow a visual interpretation ofparameter uncertainty through parameter space analysis.3.2.1 Solution to the forward problemThe finite element method was used to solve both the flow and transportmodels. For steady state flow, the finite element method is straightforward toimplement and can be solved rapidly. For transport using the standard Crank-Nicolsontime stepping scheme, numerical dispersion and stability issues constrain the size ofthe elements and the time step length required for accurate solutions. For large scaleflow systems, a large number of time steps are required for accurate solutions.40Consequently, a large amount of computer time is required in order to accuratelysimulate tracer concentrations, especially in large-scale settings The forward problemmust be solved many times during the parameter estimation process. If the computertime required to simulate tracer concentrations could be reduced, it would result in alarge overall saving in computer time required to obtain parameter estimates. In orderto be able to use parameter estimation methods in large scale flow and transportsystems, an alternate time stepping method is desirable.One alternative method for time stepping in finite element methods is theArnoldi algorithm [ Woodbury eta!., 1990]. The Arnoldi algorithm replaces the globalstiffness and storage matrices with much smaller matrices, and steps through time inthe reduced system. The finite element solution is a linear combination of Arnoldivectors, and a small number of Arnoldi vectors can capture the essentialcharacteristics of the solution. The Arnoldi vectors are calculated from the finiteelement matrices, and combined to create the reduced system of equations. Thereduced system of equations is marched through time, and a reduced solution isdetermined at each time step. At times of interest, the reduced solution is expandedinto the full solution. The numerical formulation of the Arnoldi algorithm is outlined insection 3.2.1, and the details of the Arnoldi algorithm for the transport equation aredescribed in Woodbury et a!., [1990] and Nour-Omid et al. [19911.For most problems, less than 15 Arnoldi vectors are needed to capture the41essential characteristics of the solution. Once all Arnoldi vectors are calculated, thesmall (less than 15 x 15) system is marched through time using a Crank-Nicolsonscheme. The computer time required to march this small system through time isnegligible. For a small number of time steps (less than 15) the Arnoldi method isslower than the Crank-Nicolson method. However, for many time steps, the Arnoldimethod can be much faster than the Crarik-Nicolson method.In theory, the computer time required for the calculation of a single Arnoldivector should be equivalent to a single Crank-Nicolson time step. However, thisequivalence is not true when iterative methods are used to solve the system ofequations. For the Crank-Nicolson method, after the initial time step, the starting pointfor the iterative solution at each time step is close to the true solution. The iterativemethods converge rapidly, generally in less than 5 iterations for reasonable timesteps. When calculating the Arnoldi vectors, the starting guess for the iterativemethods is far from the solution, and an average of 25-50 iterations is required forconvergence, more if the system of equations is poorly-conditioned. The algorithmsin this work use the ORTHOMIN conjugate gradient method (Mendoza et a!., 1992)to solve the systems of equations. Using ORTHOMIN, the Arnoldi method really onlystarts to outperform the Crank-N icolson method when more than 50 to 100 time stepsare required.423.2.2 Solution to the inverse problemNon-linear regression is a common method for parameter estimation. Theregression technique and solution methods used in this thesis follow Cooley and Naff[1985] and Hill [1 992]. Aset of observations Y1, I = 1,n of the physical system are fitto a model which has the parameters B1, I = 1,p. The parameters B. are the trueunknown values of the parameters. The functional relationship between theobservations and the model parameters isY=f(B)+E (3.1)where e is a vector of random variables assumed to have a distribution ofN(O, 2 V) (3.2)where V is a known (n x n) positive definite variance-covariance matrix. The matrixV defines the relative covariances between the data, and o is a scaling term with anunknown magnitude.The parameter estimation problem for joint data sets actually uses three typesof data; hydraulic head data, mass concentration data, and prior estimates of theparameters. Equations (3.1) and (3.2) can be applied to both the both the head andconcentration data, though each data set will have different functional relationshipsand error terms. The functional relationship f(B) for the head data is the steady stategroundwater flow equation (1.1), while the functional relationship f(B) for theconcentration data is the mass transport equation (1.3). For the prior information on43model parameters, the functional relationship between the observations and modelparameters isY=B-i-E (3.3)and has a distribution defined by equation (3.2). In real world problems, the use of(3!2) for the errors in the prior information is not always valid. Chapter 5 examines theconsequences of biased prior information.The non-linear regression approach minimizes an objective function S(b) toestimate the parameters b, with the resulting parameter estimates defined as 6. Forthe joint parameter estimation problem using head, concentration, and priorinformation on parameters, the objective function is:S(b)—wh(Yh—f(b))TV,1 (Yh—f(b)) + w(Y_f(b))T V1 (Y—f(b))+ (bb)T V, (b,,—b) (3.4)where the h subscript denotes those terms relating to head data, the c subscriptdenotes those terms relating to concentration data, and the p subscript denotes thoseterms relating to prior information on parameters. The Wh and w, terms are weightsdesigned to yield a trade-off between the parts of the objective function. The designand magnitude of these weighting terms are discussed in chapter 6.A modified Gauss-Newton method will be used to determine the parameterestimates for the objective function given in (3.4). The basic Gauss-Newton method44and various modifications are given in Cooley and Naff [1985]. For this work, theGauss-Newton method is modified by scaling the approximate Hessian matrix,augmenting the Hessian matrix using Marquardt’s modification [Marquardt, 1964], andsingular value decomposition (SVD) [Press et al., 1986]. This combination ofmodifications is used not only to reduce problems associated with ill-conditionedparameter sets, but also to allow easy analysis of the uncertainties associated withthe parameter estimates.The basic Gauss-Newton method is based on computing, at each iteration, theminimum of the second order Taylor expansion of the objective function (3.4). Theresulting set of equations are:Hdk=g (3.5)whereH = approximate Hessian Matrix = XT V’ X (3.6)g = Gradient Vector = XT V (Y- f(b)) (3.7)The approximate Hessian matrix (also called the least squares coefficient matrix) isa first order approximate to the true Hessian matrix. The vector dk is the correction tothe parameter vector at the kth iteration. The matrix V’ is the inverse of the combineddata reliability matrix. The matrix X is the sensitivity matrix, whose elements arewhere {f(b)}1 is the ith data point and b1 is the jth parameter. Because theobjective function (3.4) is non-linear with respect to the parameters, more than oneiteration is required to obtain the parameter estimates. A starting guess for the45parameters is used to initiate the procedure. The approximate Hessian matrix and thegradient vector are calculated using the initial estimates of the parameter values.Equation (3.5) is solved for the parameter corrections, and an updated set ofparameter estimates is calculated. These updated parameter values are then used torecalculate the approximate Hessian matrix and gradient vector, and new parametercorrections are calculated. This process is continued until the parameter valuesconverge within a given tolerance. The final parameter values are the parameterestimates. The calculation of the uncertainties in the parameter estimates, and theirrelationship to the response surfaces, are discussed in Chapter 4.3.2.2.1 Calculation of the sensitivity matrixFor the Gauss-Newton method, calculating the sensitivity matrix of partialderivatives is the most time consuming portion of the parameter estimation process.Three methods are in use [Yeh, 1986].1. Influence coefficient method. Each parameter is perturbed from its original estimateby a small amount, and a complete simulation is performed with the perturbedparameter. Since the method requires perturbing the parameters one at a time, themethod requires p+1 forward simulations to determine the sensitivity matrix.2. Sensitivity equation method. The governing equation is differentiated andrearranged so that the partial derivatives with respect to each parameter are used as46the dependent variable. The numerical solution to the set of equations yields onecolumn of the sensitivity matrix. Again, p+1 simulations are needed to determine thesensitivity matrix.3. Adioint state method. The adjoint state equations for the original partial differentialequations must be derived. The number of adjoint equations to be solved to determinethe sensitivity matrix is equal to the number of observations n.The influence coefficient method is the least accurate method [Sun and Yeh,1 990a]. The sensitivity equation method is efficient if there are more observations thanparameters to be identified. The adjoint state method is more efficient if there aremore parameters to be identified than observations. For this work, there are alwaysmore data than parameters, so the influence coefficient and sensitivity equationmethods are the best options. Both methods are used in this work, but in general thesensitivity equation method is superior, for reasons given below.Sensitivity equation method aiplied to steady state flow equationThe discretized form of the flow equation is:Ah=f (3.8)A is the (I x I) global stiffness matrix, where I is the number of nodes in the finiteelement mesh, h is the (1 x 1) vector of hydraulic head values at each node of thefinite element mesh, and t is the (I x 1) forcing vector, which includes boundary47conditions and flux terms. The sensitivity formulation of the flow equation is:(i9)abj ab1 ab1where b is the fth parameter. The variable to be solved for is , the sensitivityvector of the head to the jth parameter. The form of the equation is identical to theoriginal discretized flow equation, and only the right-hand side changes for eachparameter. The matrix consisting of the derivative of the stiffness matrix with respectto each parameter (second term on the RHS) is very sparse. Because only a sparsematrix needs to be formed for each parameter, and only the RHS of the system ofequations is altered, the sensitivity equation method is very efficient. It is much fasterthan the influence coefficient method for calculating the sensitivity coefficients in thesteady state flow equation.Sensitivity equation method applied to the transient transport equationThe discretized form of the transport equation is:[A + (1/At) M] c1 = [(1/At) M] C1 + f (3.10)A is the global stiffness matrix for the transport equation (different from the flowequation), M is the global storage matrix, At is the time step length, c is the vectorof concentrations at the new time step, c is the vector of concentrations at theprevious time step, and f is the forcing vector. The sensitivity formulation of thetransport equation is:48[A+.IM]t+At + — . {c}÷ (3.11)The variable to be solved for at each time step is , the partial derivative ofconcentration with respect to the fth parameter at the current time step. From thesolution to equation (3.10), the concentrations at each time step are known. Thederivative of the concentration with respect to the parameter isknown at the previoustime step. The form of the sensitivity equation is identical to the discretized form of thetransport equation, only the right hand side is different. However, unlike the flowequation, the matrix of the derivatives of the stiffness matrix with respect to theparameters (third term on RHS) is a full matrix. Any change in the flow parametersalters the entire velocity field, on which the global stiffness matrix is based. For thetransport equation, the sensitivity equation approach and the influence coefficientapproach are approximately equal in terms of computer time required.Sensitivity equation method applied to the Arnoldi transport equationSince the Arnoldi algorithm is used for time stepping in the transport equation,the sensitivity matrix equations need to be derived with respect to the Arnoldiequations. The Arnoldi algorithm reduces the discretized transient transport equation(3.10) to a reduced set of equations:[I + (1/At) P1 = [(1/At) P] w1 + g (3.12)where:49I = Identity matrix= [QT M A1 MO]; (mxm)g = 0T M A1 f ; (mxl)w= QTC ;(mxl)0 = matrix of Arnoldi vectors; (Ixm)m = number of Arnoldi vectorsSince the sensitivity equation is of the same form as the transport equation, thesensitivity equation (3.11) can also be reduced to the Arnoldi form. For eachparameter b., the vector and the matrix . must be formed. The vector...!.3b1 abjcan immediately be reduced into Arnoldi dimensions by pre-multiplying by (QT M A1).However, the matrix must be multiplied by the concentrations at each time step,and the resultant vector then reduced to Arnoldi dimensions and added to the RHSvector. This operation, performed at each time step, is inefficient. The advantage ofthe Arnoldi reduction is wasted.An approximate method for calculating the sensitivity coefficients using theArnoldi algorithm was developed. Instead of multiplying the matrix. by theconcentrations and then reducing the resultant vector at every time step, themultiplication and reduction is only done a few times during the entire time steppingscheme. The resulting sensitivity coefficients are not exact. The accuracy of thesensitivity coefficients depends on the number of times the multiplication and reductionis performed. One multiplication and reduction, using only the concentrations at the50end of the simulations, was generally not adequate. The convergence rate using thesesensitivity coefficients was relatively slow. Using two multiplications and reductions,the sensitivity coefficients were generally adequate. More than five multiplications andreductions resulted in good approximations to the true sensitivity coefficients for mostproblems in this thesis.Near the minimum of the objective function, especially in poorly conditionedproblems, the sensitivity coefficients need to be quite accurate. Though increasing thenumber of time steps in which the true RHS vector is calculated does increase theaccuracy of the sensitivity coefficients, it is at the expense of computer time. Thealgorithm used in this work switches from the sensitivity coefficient formulation for thetransport problem to the influence coefficient method once the convergence rate slowssignificantly. It was found that the influence coefficient method, combined with theArnoldi algorithm for transport, performs well near the minimum. Forgroundwaterflow,the sensitivity coefficient formulation is used at all iterations, since it is much moreefficient than the influence coefficient method.3.2.2.2 Modifications to the basic Gauss-Newton methodTo speed convergence for unstable problems and to allow consistent analysisof parameter uncertainties, three modifications to the basic Gauss-Newton methodwere made. The first is the Marquardt modification, which moves the descent direction51from a Newton direction toward a steepest descent direction if a reduction in theobjective function does not occur at any iteration [Marquardt, 1964]. Theimplementation of this modification is described in Cooley and Naff [1985]. The secondmodification is scaling the approximate Hessian matrix and gradient vector by theparameter values. The third modification is to use singular value decomposition (SVD)to solve the Gauss-Newton system of equations while discarding small singularvalues.The approximate Hessian matrix and the gradient vector are defined inequations (3.6) and (3.7). The elements of the approximate Hessian matrix and thegradient vector are scaled by the current parameter values.HIk = HJk * b * bk (3.13)g. g.*b. (3.14)This scaling is not the usual scaling performed in modified Gauss-Newton methods.Scaling the approximate Hessian matrix by the parameter values produces a matrixthat is independent of the differences in parameter values. Later in this work, thescaled Hessian matrix is used to characterize the relative uncertainties of theparameters as they relate to the parameter space. If the approximate Hessian matrixwere left unscaled, the parameter space would reflect not only the relative parameteruncertainties, but the differences in parameter values as well.The third modification is to use singular value decomposition (SVD) to solve for52the parameter corrections at each iteration. The scaled approximate Hessian matrixis decomposed by SVD into eigenvectors and eigenvalues (for a square matrix,elgenvalues are identical to singular values).H = U*A*UT (3.15)where U is the matrix of elgenvectors and A is the diagonal matrix containing thesingular values. The parameter corrections at each iteration are calculated by:Ab = UK1 UT * g (3.16)However, before the parameter corrections are calculated, the singular values of thescaled Hessian matrix are examined. If there are one or more singular values that aremuch smaller than the other singular values, these small singular values are set tozero. This procedure prevents the estimation routine from searching along very flatdirections of the parameter space for a minimum, and allows the routine toconcentrate on the steepest directions in parameter space. For poorly conditionedproblems, this procedure speeds convergence and allows a set of parameters to beestimated where ordinary methods would fail to converge. The singular values arealso used to characterize the shape and orientation of the parameter confidenceregions, as described in Chapter 4.53CHAPTER 4. CONCEPTS OF PARAMETER SPACE ANALYSISParameter space analysis can be a valuable tool for understanding andcalibrating groundwater flow and mass transport models. The parameter spacecontains information about the parameter estimates and their associated uncertainties,and also the shape and orientation of the confidence region. The shape andorientation of the confidence region can provide information about the interaction ofthe model parameters and the data. This chapter introduces the concepts needed forparameter space analysis. Parameter spaces and response surfaces are defined, andthe relationship between the response surface and the parameter confidence regionis formalized. The linear approximation to the parameter confidence region isexamined along with its representation by eigenvectors and elgenvalues. Examplesof response surfaces and confidence regions for well- and poorly-conditionedparameter sets are shown.4.1 Parameter spaceThe parameter space for a model is the p-dimensional space that contains allpossible combinations of parameter values, where p is the number of modelparameters. Each model parameter is represented by an axis in parameter space. Ifno information is available about the model parameters, any point in the parameterspace is equally likely as the set of parameter estimates. The process of parameter54estimation identifies the most likely region of parameter space for the estimatedparameters.In this work, the parameter spaces are scaled. Each axis in parameter spaceis scaled by a value representing that parameter, usually the parameter estimate. Thescaling is performed so that the parameter space represents the relative uncertaintiesin the parameter estimates, and not just differences in parameter values. For instance,if parameter b1 has a value of 2000, and parameter b2 has a value of 0.001, a 10%uncertainty in b1 will be much larger than a 10% uncertainty in b2. The unscaledparameter space creates the impression that b1 is much more uncertain than b2, whilethe scaled parameter space shows that b1 and b2 have similar relative uncertainties.Scaling parameter space is equivalent to comparing parameter coefficients of variationrather than comparing parameter variances.4.2 Response surfacesUsing non-linear regression, an objective function is minimized to obtainparameter estimates and their associated uncertainties. For a joint parameterestimation problem using both hydraulic head and mass concentration data, andincluding prior information on the parameters, the objective function is given byequation (3.4). Values of the objective function can be mapped into the parameterspace, and these values define a p-dimensional surface called the response surface55(also known as the objective function surface). The global minimum in this responsesurface defines the set of the parameter estimates.A response surface can be constructed for a model when the followingconditions have been met. First, a groundwater model needs to be defined andconstructed. The location and type of boundary conditions, initial conditions, anddistributed parameters need to be defined and a numerical model constructed.Second, calibration data and/or hydrogeological data need to be collected andavailable for parameter estimation. The process of parameter estimation simplysearches the response surface for a minimum.A unique response surface exists for each data value used to calibrate themodel. The shape and orientation of the response surface can indicate whatinformation the data has about the model parameters. The surfaces for each datavalue are superimposed to obtain the total response surface. The response surfacefor each data value, or each group of data values, can be analyzed. It is easy todiscover whether data values contain similar or different information about the modelparameters.When only two model parameters are estimated, maps of the responsesurfaces can be constructed as follows. A two dimensional grid is defined, with oneparameter represented by each axis. A lower and upper bound for each parameter56value and a step size for each parameter is selected. At each point on the grid (eachcombination of parameter values b1 and b2), a forward simulation is run, and thesimulated values for the hydraulic heads and tracer concentrations are calculated.These simulated data values and the parameter values are substituted into equation(3.4), and a value for the objective function is calculated. The response surface canbe visualized by contouring the values of the objective function throughout the gridarea.When more than two model parameters are estimated, it is more difficult toconstruct response surfaces. The process of constructing response surfaces withmore than two parameters makes heavy demand on computer time, since the numberof forward simulations needed increases as x, where p is the number of parameters.It is also very difficult to visualize the response surfaces in parameter spaces greaterthan two dimensions. We will present response surfaces for two parameter dimensionsonly, and use other techniques to visualize response surfaces for more than twoparameter dimensions.Using response surfaces, the estimated parameter values and associateduncertainties are easy to visualize. The Gauss-Newton method described in Chapter3 allows both the minimum in the response surface and the shape of the responsesurface to be approximated without calculating the entire response surface. Thematrices involved in the Gauss-Newton method can also be used to visualize57response surfaces for multi-parameter estimation problems.4.3 Relationship between response surfaces and confidence regionsContours of the response surface represent boundaries of the confidenceregion for the parameter estimates [Ratowsky, 1984]. Under the likelihood ratiomethod, and the assumptions that the uncertainties in the data are normallydistributed, known and included in Vh, V, and V, and that model error is negligible,an approximate (1-o)1OO% confidence region for the parameter set B is given by:S(b)-S()<2 (p) (4.1)-where S(b) is the value of the objective function using the parameter set b, definedby equation (3.4), S(B) is the value of the objective function at the estimatedparameter set B, and x2(P) is a chi-square distribution with p degrees of freedom. The? term is often unknown, and can be estimated from the model. The estimate of a2is 2, defined as:s2 S(s) (4.2)(n-p)where n is the number of data and p is the number of parameters. An approximate(1-a)100% confidence region for the parameter set B is given by:58S(b)-S() pF(p,n-p) (4.3)where F(p,n-p) is an F distribution with p and (n-p) degrees of freedom. For large n,s2 approaches a2 and pF(p,n-p) approachesx2(p), so the two confidence regions areequivalent,Figure 4.1 illustrates a response surface and confidence regions for a twoparameter system. The response surfaces are presented as (S(b)-S())Is2so that theconfidence regions can be visualized. The parameter axes b1 and b2 are scaled bytheir estimated values, so the point in parameter space at b1=1 .0, b2=1 .0 representsthe parameter estimates. Both univariate and joint confidence regions can becalculated from this response surface. A univariate confidence interval is theconfidence interval for one parameter at a time, assuming the values of the otherparameters remain fixed. For large n, the univariate confidence interval can beobtained from the 2(1) distribution. At a 68.6% confidence level (one standarddeviation), c=.314, and 2314(1)=1.0, so the univariate confidence intervals arecalculated from the 1.0 contour. The joint confidence region is the confidence regionfor both parameters, where both can vary simultaneously. For a two parametersystem, the joint confidence region is obtained from the 2(2) distribution. At a 68.6%confidence level,2314(2)23 so the joint confidence region is based on the 2.3contour. Both the univariate and joint confidence regions are shown in Figure 4.1.594.4 Linearized confidence regionsThe confidence region defined by the response surface reflects the true nonlinearity of model parameters. It is often convenient to work with a linearizedapproximation to this non-linear confidence region. The linearized confidence regioncan be calculated from the matrices involved in the Gauss-Newton parameterestimation scheme. The linearized confidence region is defined as(b_)T 5Tv12 (b-13) ps2F(p,n-.p) (4.4)where is the sensitivity matrix evaluated at the parameter estimates [Ratowsky,1984]. The boundary of (4.4) forms an ellipsoid in parameter space centered onThis ellipsoid in p-dimensional space is often hard to visualize. It is moreconvenient to calculate the major axes of this ellipsoid. The length and orientation ofthese axes can be calculated from either the sensitivity matrix, the approximateHessian matrix, or the variance-covariance matrix evaluated at the final parameterestimates. For this work, these axes are calculated from the approximate Hessianmatrix because it is readily available in scaled form, and can be manipulated byadding and subtracting data without changing dimension. The scaled Hessian matrixis decomposed by SVD into singular values and associated vectors (equivalent toeigenvalues and eigenvectors for a square matrix). The square root of the inverse ofthe eigenvalues are the lengths (L) of the axes of the ellipsoid. The associatedeigenvectors are unit vectors (U) defining the orientation of the axes of the ellipsoid60with respect to the parameter axes.The boundary of the ellipsoid is calculated to coincide with the 1.0 contour ofthe response surface. The projection of the ellipsoid onto each parameter axis definesthe univariate confidence intervals for the estimated parameters at one standarddeviation. Since the parameter space is scaled by the estimated parameter values,the univariate confidence intervals are for the coefficient of variation of the estimatedparameters.For the response surface shown in Figure 4.1, a linearized confidence regioncan be calculated from the scaled approximate Hessian matrix. The SVDdecomposition of the scaled approximate Hessian matrix yields eigenvectors andeigenvalues, which can be converted into lengths and orientation of the axesassociated with the linearized confidence ellipse. Table 4.1 summarizes the lengthsand orientation of the axes of the confidence ellipse.Table 4.1 Axes of confidence ellipsoidLengths of Axes L1 = .008 [ L2 = .024Unit Vectors U1 U2b1 axis .98 .16b2 axis .16 -.9861Constructing the ellipse from the Information contained in Table 4.1 Isstraightforward. The ellipse has two axes, a short axis with a length L1 — .008, and along axis with a length L2 .024. The unit vector U1 defines the orientation of the shortaxis with respect to the parameter axes b1 and b The unit vector U1 has a projectionof 0.98 In the direction of parameter axis b1, and a projection of 0.16 in the directionof parameter axis b2. To plot the short axis, the total projection in each parameter axisdirection is required. The total projection is the product of the axis length and the unitvector In the direction df each parameter axis. The total projection of the short axis inthe b1 axis direction is (0.008) * (0.98) —0.0078. The total projection of the short axisin the b2 axis direction is (0.008) * (0.16) — 0.001 3. The construction of the short axisusing these total projections Is shown In FIgure 4.2. For the long axis, the constructionIs similar. The total projection of the long axis in the b1 axis direction is (0.024) * (0.16)— 0.0034. The projection of the long axis in the b2 axis direction is (0.024) * (0.98) —0.023. These total projections are used to construct the long axis of the confidenceellipse, as shown in FIgure 4.2. This confidence ellipse coincides with the 1.0 contourof the confidence region.The ratio of the lengths of the longest axis of the confidence region to theshortest axis of the confidence region is termed the condition number (ON) for theparameter estimates [Sorooshlan and Gup 1985]. The ON Is a measure of thestabilIty of the parameter estimates, with smaller CNs Indicating more stableparameter estimates. For this parameter set, the ON — 3.0.624.5 Examples of response surfacesResponse surfaces are very useful in understanding and diagnosing problemsassociated with ill-conditioned parameters. The axes of the confidence regions canyield similar information. Figure 4.1 and 4.3 are examples of two parameter responsesurfaces. Figure 4.1 is a response surface for a well-conditioned parameter set, whileFigure 4.3 shows a more poorly-conditioned parameter set.For the well-conditioned set of parameters (Figure 4.1), the response surfacehas a well-defined minimum. The joint parameter confidence region enclosed by the2.3 contour is relatively small. The condition number, CN = 3.0, is relatively small. Theresponse surface for the poorly-conditioned parameter set (Figure 4.3) contains avalley with a very flat bottom. The minimum of the response surface lies somewherealong this valley. For poorly-conditioned problems, this flat-bottomed valley is thecause of small errors in data values leading to large errors in estimated parametervalues. Small errors in the data values can move the minimum in the responsesurface to anywhere along the bottom of the valley, covering a wide range ofparameter estimates. The joint confidence region for the parameter estimatesenclosed by the 2.3 contour is relatively large. The condition number for thisparameter set is large (CN = 17), showing that one axis of the confidence region ismuch longer than the other.63When more than two parameters are estimated, it is difficult to visualizeresponse surfaces. The axes of the confidence ellipsoid are used to understand theshape and orientation of the response surface. The CN can give an indication of theconditioning of the parameter set. As the parameter dimension increases, the CN willincrease, so no absolute scale can be assigned to the magnitude of the CN A:bloserinspection of the parameter axes will give more detailed information. For instance, ifall the axes of the confidence region are of similar length, then the parameter set isgenerally well-conditioned, and the CN will be small. If one axis of the confidenceregion is much longer than the other axes, a flat-bottomed valley exists in the pdimensional response surface. If two axes of the confidence region are much longerthan the other axes, a plane in p-dimensional space has values of the responsesurface near the minimum. The last two examples both represent poorly conditionedparameter sets.The orientation of the confidence ellipsoid relative to the parameter axes is alsoimportant. The orientation of the confidence region will determine how an error in theestimate of one of the parameters translates into errors in the other parameters. Thisconcept is discussed in Chapter 5. The orientation of the confidence ellipsoid is alsoimportant when using joint data sets, as the difference in orientations of theconfidence regions for the two data sets will determine the reduction in parameteruncertainty for the joint estimates. These concepts are discussed in Chapter 6.644.6 Relative contribution to the coefficient of variationTo better understand the confidence region, it is helpful to know the contributionfrom each axis of the confidence ellipsoid to the coefficient of variation of theparameter estimates. The total coefficient of variation for each parameter, CV, is thesquare root of the element (i,i) of the inverse scaled approximate Hessian matrix. CV1can also be calculated from the axes of the confidence ellipsoid asCV=lJj=1,p (U1L)2 (4.5)Conceptually, the squared projection of all axes of the confidence region j in thedirection of the parameter i are summed. The total coefficient of variation is the squareroot of this sum. Figure 4.4 is the confidence ellipse for the response surface of Figure4.1. The axes of the confidence region and the contribution of each axis to the totalcoefficient of variation of parameter b1 are labeled. The relative contribution to the totalcoefficient of variation of each parameter i from each parameter axis j is‘U L’RC.. = ‘ ‘‘ “ (4.6)CThe relative contribution (RC) has values ranging from zero to one. If a relativecontribution RC is near zero, then ellipsoid axis j has a very small contribution to thetotal coefficient of variation of parameter i. If a relative contribution RC is near one,then ellipsoid axis j has a very large contribution to the total coefficient of variation ofparameter I. The relative contributions for the axes of the confidence ellipse in Table654.1 are shown in Table 4.2.Table 4.2 Relative contribution to total CV from each axis of the confidence regionAxis 1 Axis 2Parameter b1 .81 .19Parameter b2 .01 .994.7 SummaryThis chapter introduced the concepts needed for parameter space analysis.Response surfaces are a picture of the objective function in parameter space. Bothunivariate and joint confidence region for the parameter estimates can be determinedfrom the response surface. A linearized approximation to the confidence region maybe calculated from the approximate Hessian matrix, and SVD used to calculate theaxes of the confidence ellipsoid. The total coefficient of variation and the relativecontribution to the coefficient of variation from each axis of the joint confidenceellipsoid can be determined. Examples have been presented for well- and poorlyconditioned parameter sets.66Univariate confidenceinterval for b2Joint confidenceregir for b1— L2Figure 4.1 Response surface, univariate confidence intervals and joint confidenceregion for parameters b1 and b2b10.95 1.0 1.05N1.05b2 1.00.95N N18.4Univariate confidence inteaI for b167b2 1.0.0078Axis 1‘ 0013.023\)jAxis 2Figure 4.2 Construction of the axes of the confidence ellipse from the axis lengthsand unit vectors0.951.0501 1.050.95-‘.0034680.95 1.051.05b2 1.00.95Figure 4.3 Example of a response surface for a poorly-conditioned parameter set.69b2 1.0Figure 4.4 Axes of confidence region showing partial coefficient of variation fromeach axis of the confidence region.700.951.05b1•L0 1.050.95CHAPTER 5. PARAMETER SPACE ANALYSIS FOR EFFICIENT ANDRESPONSIBLE USE OF PRIOR INFORMATIONParameter space analysis can be used to develop guidelines for efficiently andresponsibly stabilizing the inverse problem using prior information. When a flow modelis constructed, all parameters associated with the model need to be defined, both theinternal model parameters such as hydraulic conductivity and recharge, and theboundary values such as specified head and specified flux. If all model parametersare estimated using only hydraulic head data, the inverse problem will be non-identifiable. This non-identifiability can be observed in the response surface for theparameter set. The response surface will contain a perfectly flat-bottomed valleycutting across the entire parameter space.To change this non-identifiable inverse problem into an identifiable problem,prior information is needed. This prior information can take one of two forms: (1) aspecified value for a given parameter, assuming no uncertainty in the specified value,or (2) a prior value with some uncertainty. Specifying a parameter value with nouncertainty is equivalent to eliminating a dimension in parameter space. Priorinformation with uncertainty alters the topology of the response surface in the fullparameter space. These two uses of prior information are treated separately in thischapter because they act differently in parameter estimation problems and producedifferent results in parameter space. In this work, prior information without uncertainty71is referred to as specifying a parameter value. Prior information with uncertainty istermed prior information.In many cases, when either of the above approaches is used to eliminate thenon-identifiability of the parameter set, the parameter set remains unstable. Thisinstability can be observed in the response surface for the parameter set. Theresponse surface will contain a valley with a nearly flat boffom. The unstableparameter set is identifiable, since a unique minimum does exist, but it is very difficultto find this minimum. The response surfaces for non-identifiable and unstableparameter sets are very similar. The same remedies used to correct problems withnon-identifiable parameter sets can be used to correct problems with unstableparameter sets.When prior information is used to stabilize a set of parameter estimates, it mustbe recognized that the prior information may be in error with respect to the model. Asdiscussed in Chapter 2, the prior information may be biased. If errors in the priorestimates are present, the final parameter estimates for all parameters may beinfluenced by the erroneous prior information. The influence of errors in the priorinformation on the parameter estimates can be determined using response surfacesand parameter space analysis.In this chapter, guidelines will be developed to ensure the most efficient and72responsible use of prior information. Those parameters for which prior information willreduce the overall uncertainty in the parameter estimates the most and produce themost stable parameter set will be the most efficient parameters. Responsible use ofprior information involves identifying how errors in the prior information will influencethe final parameter estimates. Those parameters for which errors in the priorinformation will have the least influence on the final parameter estimates will be themost responsible parameters. Both the most efficient and most responsibleparameters can be identified through an examination of the model parameter spaceand response surfaces.This chapter is organized as follows. First, a synthetic flow model is introduced.The parameter space of this synthetic model is used to illustrate the conceptsdeveloped in this chapter. The effect of specifying parameter values is explored usinga two parameter subset of the synthetic model. The effect of prior information isexplored using the same two parameter subset. Guidelines for identifying thoseparameters for which prior information is most efficient are developed. Next, theconsequences of error in the prior parameter values is explored using a threeparameter subset of the synthetic model. Guidelines for identifying those parameterswhich result in the most responsible parameter estimates are developed. Finally, theconcepts and methods developed using the simple two and three parameter systemsare demonstrated on a more complicated multi-parameter system.735.1 Synthetic modelA simple synthetic flow model is used to illustrate the concepts developed inthis chapter. A two-dimensional flow system is constructed. The flow system is 6 kmby 6 km, and contains two transmissMty zones, one zone of enhanced areal recharge,and two constant head boundary zones (Figure 5.1). The remaining model boundariesare no-flow boundaries. The flow system is at steady state. Fifteen observation pointsare used. Observed values of the hydraulic head are simulated by running a forwardsimulation using the true parameter values, and calculating the data values at thesampling locations. Random Gaussian errors with an uncertainty of 1 meter are addedto these data values to generate the observed data values. These observed datavalues are used to estimate the parameters.There are five parameters in the flow system: the two transmissivity zones (T1and T2); one recharge zone (R); and two specified head zones (H1 and H2). Duringparameter estimation, the transmissivity parameters are log-transformed, and the otherparameters are untransformed. For ease of presenting response surfaces, generallyonly two parameters will be estimated at one time. The other three parameters arefixed at their true values. In the following sections, the effects of including priorinformation are examined using the parameter space for parameters T1 and R. Thisparameter space is examined using both the response surface and the axes of theconfidence region.745.1.1 Parameter space T-RFigure 5.2 illustrates the response surface for the parameters T1 and R. Theresponse surface is plotted as [S(b)-S(6)]/s2,so the contours represent the height ofthe response surface above the minimum. Since the parameter space is scaled by theparameter estimates, the parameter estimates are at the point T1 = 1.0, R = 1.0 inparameter space. The response surface contains a long, narrow valley. The axis ofthis valley is oriented such that is makes a much smaller angle with respect to the Rparameter axis than the T1 parameter axis.The confidence ellipse is a linearized picture of the response surface, with theboundary of the ellipse corresponding to the 1.0 contour of the response surface. Thisellipse can be characterized by its major axes, calculated from SVD decomposition ofthe scaled Hessian matrix at the parameter estimates. Table 5.1 gives the axes of theconfidence ellipse for this parameter set. L1 and L2 are the lengths of the two axes,and U1 and U2 are the unit vectors showing the orientation of the axes. The conditionnumber for this parameter set is 16.8, indicating that axis 2 is approximately 17 timeslonger than axis 1.Both the response surface and the axes of the confidence ellipse indicate thatalong the direction of axis 2 the response surface is relatively flat. The parameter setT1-R is relatively unstable. Using equations (4.5) and (4.6), the relative contributions75Table 5.1 Axes of confidence ellipse for T1-R parameter setLengths of Axes L1 = .0077 L2 = .129 1Orientation of Axes U1 U2T1 axis -.94 -.32R axis .32 .94to the parameter coefficient of variation from each axis of the confidence region canbe calculated (Table 5.2). Axis 2 is the major contributor to the uncertainty to bothparameters. If the length of axis 2 could be reduced, then the uncertainty in bothparameter estimates would be decreased.Table 5.2 Relative contributions to the parameter CV from each axis of theconfidence region for parameter set T1-RAxis 1 Axis 2Parameter T1 .029 .971Parameter R .001 .9995.2 Efficient use of prior informationIn order to stabilize this parameter set, prior information can be used. The priorinformation may either be prior information without uncertainty (equivalent to specifyinga parameter value), or prior information with uncertainty. First, consider specifying aparameter value, which results in removing a parameter dimension from parameter76space. The response surface for this reduced parameter set can be visualized bytaking a cross section through the two-dimensional response surface (Figure 5.2)perpendicular to the specified parameter axis, intersecting that axis at a pointcorresponding to the specified value. The intersection of this cross section with the 1.0contour is the univariate confidence interval at one standard deviation. If parameterR is specified at a scaled value of 1.0 and parameter T1 estimated, the 68%confidence interval for T1 is approximately 0.99<T1<1.01. If parameter T1 is specifiedat a scaled value of 1 .0 and parameter R is estimated, then the confidence intervalfor R is approximately 0.97<R<1 .03. Specifying R results in a smaller confidenceinterval than specifying T1. For the purpose of reducing parameter uncertainty, it ismore efficient to specify R than to specify T1. This simple example demonstrates thatis it most efficient to specify the parameter whose axis is most closely aligned with thevalley in the response surface.When one parameter is specified, and thereby removed from the inverseproblem, that parameter is assigned a single value with no uncertainty. However, theparameter whose axis is most closely aligned with the valley in the response surfaceis often the one with the most uncertainty associated with it. The above analysissuggests that using prior information on this parameter will most efficiently stabilizethe remaining parameter estimates. This can be puzzling; the parameter with the mostuncertainty is replaced with a parameter with no uncertainty. It is the structure of theresponse surface which makes this demand- based on the data available, the77response surface contains the least information about this parameter. If thisparameter, about which the response surface space contains very little information,does not need to be estimated, the other parameters can be estimated moreprecisely. However, there is an obvious danger in specifying a value of a parameter.If the value of the specified parameter is in error/the estimate of the remainingparameter will also be in error, and this error Will not be accounted for in thecalculation of the confidence region. This issue is discussed in section 5.4.Prior information with uncertainty can also be used to stabilize the parameterset. Figure 5.3 illustrates the response surface for prior information on parameter Ronly, where the scaled prior value equals 1 .0 with an uncertainty of 2%. The totalresponse surface for head and prior information is shown in Figure 5.4. The size ofthe confidence region is reduced, and the parameter set is more stable than usinghead data only. To demonstrate the efficiency of prior information on differentparameters, the lengths of the longest axis of the confidence region after includingprior information can be compared. The length of the longest axis of the confidenceregion using only head data was 0.129 (from Table 5.1). For head and priorinformation on T1, the length of the longest axis is 0.068. For head and priorinformation on R, the length of the longest axis is 0.031. Prior information on R ismost efficient, since it reduces the length of the longest axis of the confidence regionthe most with a given quality of prior information.78These graphical approaches using response surfaces can be formalized usingthe axes of the confidence region. Prior information on the parameter with the largestelement of the longest axis of the confidence region will reduce the uncertainty of theremaining parameter estimates most efficiently. The parameter with the largest- - element of the longest axis of the confidence region is the parameter whose axis ismost closely aligned with the longest axis of the confidence region. From Table 5.1,axis 2 is the longest axis, and parameter R has the largest element of that axis (0.94).Prior information is most efficient for those parameters with the largest elements of thelongest axis of the confidence region.5.3 Effect of errors in prior informationUsing prior information is an effective way to stabilize an unstable or nonidentifiabe inverse problem. However, if errors exist in the prior information, theseerrors will influence the estimates of the remaining parameters. Prior informationwithout uncertainty is used to illustrate the consequences of errors in prior informationand develop the methods of minimizing the consequences of these errors. Thesemethods are then extended to using prior information with uncertainty.Based on examining a number of cases, it is apparent that in some cases anerror in the specified value of a model parameter leads to very large errors in theestimated parameter values of the other model parameters. In other cases, an error79in the specified value of a model parameter value leads to relatively small errors in theestimated parameter values of the other model parameters. The ratio of a scaled errorin an estimated parameter value to a scaled error in a specified parameter value istermed the error ratio. This section examines the causes of the error ratio, and thereason for different error ratios for different model parameters. The underlying factorthat determines the magnitude of the error ratio is shown to be the orientation of theconfidence region in parameter space, which leads to a method of approximating theerror ratios.5.3.1 Consequences of errors in prior estimatesTo illustrate the consequences of error in prior estimates, the flow systemintroduced above is used, and parameters T1 and T2 are estimated using hydraulichead data. The data has an uncertainty with a standard deviation of 1 meter. In thissystem, there are three other parameters that have been assigned specified values.These are the inflow head boundary, the outflow head boundary, and the rechargezone. When the parameters T1 and T2 are estimated using only the head data, andthe true values of the other parameters are used, the true values of parameters T1and T2 are estimated. Figure 5.5 is the response surface for the parameter set T1-T2using only head data, with the confidence region due to uncertainty in the head data.If an error exists in the specified value of one of the other parameters, the80estimates of parameters T1 and T2 will not be equal to their true values. Table 5.3illustrates the consequences of an error in the inf low head boundary (H) and Table 5.4illustrates the consequences of an error in the recharge (R) parameter. For theparameter H, relatively small errors in the specified value lead to large errors in theestimate of parameter T2. The error ratio is calculated as a ratio of the percentageerror in the estimated parameter to a percentage error in the specified value.The errorratios for parameter T2 are much larger than one. The errors in the estimate of T1 areapproximately of the same magnitude as the errors in the specified head value, withthe error ratios of approximately one. When specifying R, errors in the specified valuelead to much smaller errors in the estimate of both T1 and T2. The error ratios forspecifying R are much smaller than one.Table 5.3 Consequences of error in specified head boundary valueValue of % Error Estimate % Error Error Estimate % Error ErrorH in H of T, in T1 Ratio of T2 in T2 Ratio200.00 0.00 1.290 0.00 0.00 2.231 0.00 0.00201.00 0.50 1.294 0:31 0.62 2.308 3.45 6.90205.00 2.50 1.314 1.86 0.74 2.884 29.27 11.71210.00 5.00 1.360 5.43 1.09 6.701 200.36 40.07199.00 -0.50 1.282 -0.62 1.24 2.163 -3.05 6.10195.00 -2.50 1.262 -2.17 0.87 1.958 -12.24 4.89190.00 -5.00 1.241 -3.80 0.76 1.779 -20.26 4.05180.00 -10.00 1.212 -6.34 0.63 1.531 -33.67 3.37150.00 -25.00 1.124 -12.87 0.51 1.100 -50.69 2.0381Table 5.4 Consequences of error in specified value of rechargeValue of %Error Estimate % Error Error Estimate % Error ErrorR in R of T1 in T1 Ratio of T2 in T2 Ratio0.00040 0.00 1.290 0.00 0.00 2.231 0.00 0.000.00041 2.50 1.299 0.70 0.28 2.242 0.49 0.200.00042 5.00 1.310 1.55 0.31 2.250 0.85 0.170.00044 10.00 1.330 3.10 0.31 2.273 1.88 0.190.00045 12.50 1.340 3.88 0.31 2.283 2.33 0.190.00050 25.00 1.386 7.44 0.30 2.329 4.39 0.180.00060 50.00 1.465 13.57 0.27 2.408 7.93 0.160.00080 100.00 1.590 23.26 0.23 2.534 13.58 0.140.00039 -2.50 1.279 -0.85 0.34 2.220 -0.49 0.200.00038 -5.00 1.270 -1.55 0.31 2.211 -0.90 0.180.00036 -10.00 1.245 -3.49 0.35 2.188 -1.93 0.190.00035 -12.50 1.230 -4.65 0.37 2.173 -2.60 0.210.00030 -25.00 1.173 -9.07 0.36 2.104 -5.70 0.230.00020 -50.00 0.989 -23.33 0.47 1.930 -13.49 0.270.00010 -75.00 0.688 -46.67 0.62 1.628 -27.03 0.36The effect of errors in the specified parameter values is summarized graphicallyin Figure 5.6. The parameter estimates produced when errors are present in thespecified parameter values are plotted in the scaled parameter space T1-T2. For errorsin H, the parameter estimates plot along a slightly curved line in parameter space. Theerrors are much larger if H is overestimated than if H is underestimated, and even a5% overestimate in H produced parameter estimates which plotted outside of the82diagram. For errors in R, the parameter estimates plot along a relatively straight linein parameter space, and the error ratios are nearly constant. The orientation of theline in parameter space reflecting errors in R is different from the orientation of the linereflecting errors in H.For this synthetic example, several observations can be made about theconsequences of errors in specified parameter values on parameter estimates.1. Errors in the values of specified parameters lead to errors in estimatedparameter values.2. The estimated parameter values produced when errors are present in thespecified parameter values plot in relatively straight or slightly curved lines inparameter space.3. The orientation of the line in parameter space connecting the estimatedparameter values produced when errors are present in specified parametervalues bears no relation to the orientation of the parameter confidence region.4. Errors in the specified value of recharge lead to much smaller error ratios forthe parameters T1 and T2 than errors in the specified value along the inflowhead boundary.5. The orientation of the line in parameter space based on errors in R is differentfrom the orientation of the line based on errors in H.835.3.2 Conceptual relationship between confidence region and error ratioIn order to predict the consequences of errors in prior information, theconfidence region for both the estimated parameters and the parameters with priorinformation must be examined. For example, the confidence region for T1, T2, and Hmust be examined to determine the consequences of errors in H on the estimates ofT1 and T2. A three dimensional parameter space is required. It is possible to use theaxes of the linearized confidence ellipsoid to represent the three dimensionalconfidence region.The confidence region for the parameters T1, T2, and H is examined first. SVDdecomposition of the scaled Hessian matrix using the true parameter values is usedto calculate the principal axes of the confidence region (Table 5.5). Axis 3 is thelongest axis, about five times as long as axis 2 and 50 times as long as axis 1, andU3 gives the orientation of the unit vector associated with the longest axis. Figure 5.7shows the orientation of the longest axis of the ellipsoid in the three dimensionalparameter space, with the shaded box included to show the orientation of the axis.The longest axis of the parameter confidence ellipsoid represents the trace ofthe smallest value of the response surface in three dimensional space. If all threeparameters are considered, then the overall minimum of the response surface is atthe center of the ellipsoid. To determine the parameter estimates and confidence84Table 5.5 Axes of confidence region forT1-T2H parameter setAxis Lengths 0.005 0.020 0.100Unit Vectors U1 U2 U3T1 -O.i7 0.98 -0.06T2 -0.10 -0.07 -0.99H 0.98 0.17 -0.11region for T1 and T2 at any value of H, a slice parallel to the T1-T2 plane is taken,intersecting the H axis at the specified value of H. The point where the minimum inthe response surface intersects this slice is approximately the parameter estimates forT1 and T2. This minimum in the response surface is closely related to the point wherethe longest axis of the confidence ellipsoid intersects the slice.In this example, the confidence ellipsoid is oriented in such a way that theangle between its longest axis and the parameter axis T2 is relatively small (seeFigure 5.7). The angles between the longest axis of the parameter confidence ellipsoidand the other two parameter axes are relatively large. If the value of H is equal to thetrue value, then the estimates for T1 and T2 will be the equal to their true values aswell. If the specified value of H is larger than the true value, the parameter estimatefor T2 changes significantly. Since the longest axis is aligned closely with theorientation of the T2 axis, small errors in the specified value of H result in large errorsin the value of T2. However, errors in the specified value of H are approximately equalto the errors in the estimated values of T1. Comparing Figure 5.7 to Figure 5.6, the85line of parameter estimates at different errors in H is approximately the trace of thelongest axis of the parameter confidence region in three dimensional space projectedonto the T1-T2 plane. The orientation of this axis explains the fact that small errors inH resulted in large errors in the estimate of T2.The same analysis can be carried out for errors in the recharge parameter. Thelongest axis of the parameter confidence ellipsoid in the parameter spaceT1-T2R isshown in Figure 5.8, and all axes are listed in Table 5.6. In this case, the longest axisof the confidence ellipsoid is 375 times as long as the next longest axis. The threeparameters are very unstable when considered together, and the system is ill-conditioned. The longest axis of the confidence region is oriented at a small angle tothe parameter axis R, and at large angles to the parameter axes T1 and T2. Theparameter estimates of T1 and T2 at any value of recharge can be determined by thepoint of intersection of the long axis with a slice taken parallel to the T1-T2 plane atthat value of R. The longest axis is oriented at a small angle to the axis of parameterR, so large errors in recharge result in relatively smaller errors in T1 and T2. The lineof parameter estimates in Figure 5.6 is the projection of the long axis in threedimensional space on the T1-T2 plane.The error in the estimated parameters as a result of an error in the specifiedvalue of a parameter depends on the orientation of the joint confidence region. If thelongest axis of the confidence region is nearly parallel to the specified parameter axis,86Table 5.6 Axes of confidence region forT1-T2R parameter setAxis Lengths 0.016 18.100 0.048Unit Vectors U1 U2 U3T1 -0.93 0.30 0.21T2 -0.16 0.17 -0.97R 0.33 0.94 0.11then the error ratio is small. If the longest axis is oriented nearly perpendicular to thespecified parameter axis, then the error ratio for at least one of the estimatedparameters is large.5.3.3 Calculation of error ratios from confidence regionIt is possible to calculate the error ratio from the parameter confidence region.Figure 5.9 is a generic confidence ellipse in two parameter dimensions with theparameter axes scaled by the parameter estimates. Using linearized error analysis,the shape and orientation of the confidence ellipse is constant for any contour of theresponse surface, so the confidence ellipse in Figure 5.9 can represent any contourof the response surface. The error ratio is defined as the ratio of an error in theestimated parameter value to an error in the specified parameter value. If b1 is thespecified parameter, the error ratio can be graphically constructed by using a lineperpendicular to the b1 axis tangent to the confidence ellipse. This line is representedby A - A’. The parameter estimates, given an error in b1, are at the point where the87line A - A’ intersects the confidence ellipse, labelled as point (X,Y). The objectivefunction has a minimum along the line A - A’ at that point. The error in the specifiedparameter b1 is given by (X - X0). The error in the estimated parameter is given by (Y- Y0). The error ratio is therefore (Y-Y0)/(X-X. Because the shape and orientation ofthe confidence ellipse is constant for any contour of the response surface, thedistance (X-X0)can be scaled to any error in the specified parameter. The error ratioremains constant for any error in the specified parameter.The conceptual method presented above can be used in simple two and threedimensional parameter spaces, but not for multi-parameter problems. The error ratiocan be calculated from the axes of the parameter confidence ellipse. The error ratioisER - k=l,p1.kULk’ (5.1)ii- cv12where ER11 is the error ratio for an error in the ith estimated parameter, given an errorin the jth specified parameter. CV is the total coefficient of variation for the specifiedparameter, defined in equation (4.5). This error ratio is termed the linearized error ratiobecause it is calculated from the linearized confidence ellipsoid.The error ratio depends on two factors; the orientation of the confidence regionand the shape of the confidence region. If the confidence region has a large conditionnumber (long and narrow), the error ratio will depend mainly on the orientation of thelongest axis of the confidence region. Figure 5.1 Oa shows an example of a confidence88region with a large condition number, along with the longest axis of the confidenceregion. The line perpendicular to b1 axis, tangent to the confidence region, intersectsthe confidence region near the longest axis of the confidence ellipsoid. The ratio ofthe elements of the eigenvector of the longest axis of the confidence regionapproximates the linearized error ratio. If the confidence region is well-conditioned, thecontribution to the coefficient of variation from the longest axis for the fixed parameterwill be less than one, and the error ratio will be reduced from that calculated usingonly the elements of the elgenvector. Figure 5.lOb shows an example of a wellconditioned confidence region, and the longest axis of the confidence region. The lineperpendicular to the P1 axis, tangent to the confidence region, intersects theconfidence region at a point somewhat below the longest axis of the confidenceregion. Thus the linearized error ratio is less than that calculated only using thelongest axis of the confidence region.5.3.4 Using error ratios to identify responsible parameters for prior informationWhen prior information is used to reduce the uncertainty in parameter estimatesand stabilize parameter sets, it should be recognized that errors in the priorinformation may lead to errors in the estimated parameters. It is necessary to knowthe influence of errors in the prior information on the parameter estimates in order touse the prior information in a responsible manner. If errors exist in the priorinformation, it would be responsible use prior information which minimizes the effects89of these errors. The parameters that lead to the smallest error ratios are the mostresponsible parameters for prior information.Tables 5.7 and 5.8 show the linearized error ratios calculated from equation(5.1) for the two parameter sets described above. These tables can be used tochoose the responsible parameters for prior information. For the parameter set T1-T2H, the largest linearized error ratio results from estimating T2 while parameter H isspecified. Errors in specified values of H will be magnified over 7 times when T2 isestimated. At the same time, errors in the estimate of T1 will be smaller than the errorin the specified value of H. From examining Table 5.7, it is apparent that specifyingthe value of T2 leads to the smallest error ratios for the estimates of the remainingparameters. The errors in estimates of both T1 and H will much smaller than errors inthe specified value of T2. The choice of T2 as the most responsible parameter tospecify is consistent with the examination of the parameter confidence region. Thelongest axis of the confidence region is sub-parallel to the T2 parameter axis, so errorsin T2 will not result in large error ratios.90Table 5.7 Error ratios for the T1-T2H parameter setFixed T1 Fixed T2 Fixed HEstimated T1 1.0 .057 .86EstimatedT2 1.35 1.0 7.11Estimated H .30 .111 1.0For the parameter set T1-T2R, the error ratios are shown in Table 5.8. Thesmallest error ratios result from specifying parameter R. This result is consistent withthe fact that the confidence region for this parameter set is sub-parallel to theparameter axis R.Table 5.8 Error ratios for theT1-T2R parameter setFixed T1 Fixed T2 Fixed REstimated T1 1 .0 1.7 .32Estimated T2 .57 1.0 .19Estimated R 3.1 5.4 1.0It is important to note that the error ratios for the parameter subset T1-T2 arenot always the same. From Table 5.7, an error in the specified value of T1 leads to a1.2 error ratio for parameter T2. However, from Table 5.8, an error in the specifiedvalue of T1 lead to a .57 error ratio for parameter T2. The full confidence region for theparameter set determines the error ratio. The confidence regions for these twoparameter sets are different, so error ratios for the same parameters are different.915.3.5 Influence of non-linear confidence regionsThroughout the above analysis, linear approximations to the confidence regionshave been used to calculate the axes of the confidence ellipsoids and the linearizederror ratios. These linear approximations do not accàunt for the non-linearity of thetrue confidence regions. Three indications of no-linearity are apparent in theexamples presented, especially in theT1-T2H system.1. The lines in parameter space reflecting the trace of the axis of the confidenceregion are slightly curved for the confidence region including parameter H. Thiscurvature implies that the axis of the confidence region is curved in parameterspace. For a linear model, the axis of the confidence region would be straight.2. Errors in specified parameters that overestimate the true parameter values resultin different errors in the estimated parameters than errors in specifiedparameters that underestimate the true parameter values. The confidenceregion is not symmetrical about the true parameter value. In a linear system,the confidence region would be symmetrical.3. The true error ratios are not constant as the error in specified parameter valuesincreases. A combination of curved and non-symmetrical confidence regionsare the cause of the non-constant error ratios. For a linear system, the errorratio would be constant for any error in the specified parameter value. Thelinearized error ratios are considered constant.92Even though there are indications of model non-linearity, the non-linearity isgenerally not severe. The curvature of the confidence region is slight, while the non-symmetry is somewhat greater. However, since the linearized error ratio analysisdepends mainly on the orientation of the parameter confidence region, and not on thesymmetry of this region, the linear approximation is generally robust. The true errorratios may differ from the linearized error ratios, but the conclusions based on thelinearized error ratios are valid. The parameters that result in the smallest error ratioswill most likely be the same even if the true non-linearity were taken into account.5.4 Error ratios for prior information with uncertaintyUsing prior information with uncertainty is generally an improvement over simplyspecifying a single value for a parameter. In parameter space, specifying a singlevalue is equivalent to taking a slice in parameter space perpendicular to thatparameter axis at the specified value. The ill-conditioning is reduced by reducing theparameter dimension. Prior information retains all the parameters in the inverseproblem, but adds a penalty term to the response surface for deviations from the priorestimate. A beffer conditioned, more stable response surface is the result.The error ratio method developed above can be extended to using priorinformation with uncertainty. In terms of error ratio, the major difference between usingprior information with uncertainty and specifying a parameter value without uncertainty93lies in how close the final estimate is to the prior estimate. When a parameter valueis specified, it cannot change during parameter estimation, and the final estimate isequal to the prior information. When using prior information with uncertainty, the priorestimate and the final estimate are different.Two factors determine how close the final estimate are to the prior estimate;the assigned uncertainty in the prior information and the shape of the responsesurface before prior information is included. If the uncertainty assigned to the priorestimate is small, there is a large penalty in moving away from this prior estimate, andthe final estimate will be close to the prior estimate. Conversely, if the uncertainty inthe prior estimate is large, the penalty for moving away from the prior estimate issmall.If the response surface using only head data reflects a poorly-conditionedproblem, the final estimate will be close to the prior estimate. A poorly conditionedparameter set is characterized by a valley with a flat bottom. The prior information,even if it has a large uncertainty, will produce a definite minimum in the valley at thevalue of prior estimate. The final estimate will be very close to the prior estimate. Ifthe response surface using only head data is well conditioned, the minimum definedby the head data will have a large influence on the final parameter estimates. Evenif the prior information has large errors, the final estimates will be controlled by theresponse surface for the head data.94When prior information is assigned for more than one parameter, therelationships between the prior and final estimate are more complicated. In additionto the factors described above, the final estimates also depend on the relationshipbetween the orientation of the data response surface and the location of the minimumand shape of the response surface using prior information. These relationships are toocomplicated to form any specific conclusions. In general, we find that for an inverseproblem that is ill-conditioned before adding prior information, one or more of the finalestimates will be very constrained by the prior estimate. If the final estimate is in error,it will introduce errors in the final parameter estimates of the other parameters.Thus, for an ill-conditioned or non-identifiable inverse problem, the finalparameter estimates will be close to the prior estimates, regardless of the assigneduncertainty in the prior estimates. The error ratio method will be valid for determiningthe effects of errors in the prior information. For a well-conditioned problem, the finalestimates may not be close to the prior estimates, and the error ratios may be lessthan those calculated by the linearized error ratios.5.4.1 Example of the influence of the topology of the data response surfaceResponse surfaces T1-T2 and T1-R are used to illustrate the influence of thetopology of the response surface on the difference between the prior and the finalestimate. The response surface for T1-T2 (Figure 5.5) reflects a well-conditioned95parameter set, and the response surface for T1-R (Figure 5.2) is poorly-conditioned.Assume that the prior estimate has a value of T1 (prior) = 1.02 in scaled space, andthus is in error. Figure 5.11 shows the response surface for the prior information onT1, with an uncertainty in the prior information equal to one percent of the priorparameter estimate. The final parameter estimates for the two parameter sets arequite different. The total response surfaces are shown in Figure 5.12 and 5.13. Notethat the two response surfaces are not shown at the same scale. For the parameterset T1-R (Figure 5.12), the final parameter estimate for T1 is the same as the priorparameter estimate. The response surface using head data alone does not have aunique minimum, but instead contains a nearly flat bottomed valley. All along thisvalley, the values of the response surface are very close to the minimum. The priorvalue of the parameter T1 defines where the minimum will lie along this valley. Sincethe prior estimate for T1 was in error, the final estimates for both T1 and R are in error.For the parameter set T1-T2 (Figure 5.13), the parameter estimate forT1 moves awayfrom the prior estimate toward the parameter estimate based on the head data only.The head data alone produce a response surface with a definite minimum. Thisminimum has a large influence on the total response surface, and pulls the finalestimates toward this minimum.In parameter set T1-T2, the parameter space is well conditioned, and theminimum is well defined using only head data without prior information. When priorinformation is included, the parameter estimates are pulled toward the minimum using96head data only. Even when the prior information has large errors, the final estimatesare controlled by the information contained in the head data. In parameter set T1-R,the parameter space is poorly conditioned using only head data, as shown by a long,flat bottomed valley in the response surface. The prior information does help transformthis problem into a well-conditioned inverse problem, with a well-defined minimum.However, the location of the minimum (and therefore the parameter estimates) iscompletely dependent on the specified value of the prior information. If the priorinformation is accurate, there is no problem. However, if the prior information is inerror, even if a large uncertainty is assigned to the prior information, the value of theparameter estimates depends on the incorrect value of the prior information.5.5 Summary of guidelines for use of prior informationParameter space analysis has been used to develop guidelines for the efficientand responsible use of prior information in groundwater flow models. The axes of theconfidence region can be used to identify the most efficient parameters for priorinformation. Those parameters with the largest elements of the longest axis of theconfidence region will be the most efficient parameters for reducing parameteruncertainty and increasing parameter stability.A linearized error ratio matrix can be calculated from the axes of the confidenceregion, and used to identify the most responsible prior information. The linearized error97ratio matrix is calculated based on prior information without uncertainty, but is alsovalid for prior information with uncertainty. For ill-conditioned or non-identifiable inverseproblems, the error ratio method will be valid for determining the effects of errors inprior information with uncertainty. For a well-conditioned problem, the error ratios maybe less than those calculated by the linearized error ratios.To minimize the effects of possible errors in the prior estimates, thoseparameters with the smallest error ratios are the most responsible parameters for priorinformation. Errors in these parameters will have the least influence on the parameterestimates for the remaining parameters. Those parameters that are the most efficientare often the same ones as those that are the most responsible. Both guidelines leadto selecting parameters whose axes are closely aligned with the longest axis of theparameter confidence region.5.6 Application to a multi-parameter systemThe guidelines developed above for the efficient and responsible use of priorinformation to stabilize an inverse problem were applied to a synthetic multiparametersystem. The flow system is taken from Carrera and Neuman (1986c) and shown inFigure 5.14. The data consist of 18 head measurements with uncertainties of 1 meter.The parameters to be estimated include 9 transmissivity zones, 2 recharge zones, onespecified flux boundary and one specified head boundary. All model parameters,98including those describing the domain boundaries, were estimated, which resulted ina non-identifiable parameter estimation problem. When estimating parameters, thefinal estimates were found to depend on the initial values of the parameters.Table 5.9 gives the parameter estimates and uncertainties for two sets of initialparameter values. The transmissivity parameters are log transformed, all otherparameters are untransformed. One set of initial estimates is close to the true valuesof the parameters, the other set of initial estimates is far from the true values of theparameters. The two final sets of parameter estimates are very different. The standarderrors of the estimates from the two sets of initial estimates are very similar, so onlyone set is given. The standard errors of the estimates are very large compared to thevalues of the estimates, leading to very large coefficients of variation for nearly allparameters. The correlation matrix, which is not presented here, shows correlationsof greater than .99 between all parameters except H1, which showed very littlecorrelation with the other parameters. The resulting parameter estimates are verypoor, due to the ill-conditioning of the problem.5.6.1 Parameter space of multi-parameter problemBefore using prior information to stabilize this problem, the shape andorientation of the parameter confidence region must be understood. The axes of theparameter confidence region are calculated using SVD decomposition of the scaled99Table 5.9 Parameter estimates and uncertainties for multiparameter problemParameter EstimatesParameters True Value Run #1 RUn #2 Standard CVErrorT1 2.17 2.22 1.81 5.76 3.16T2 2.17 2.09 1.68 5.66 3.36T3 2.17 2.00 ‘1.60 5.31 3.32T4 2.17 2.28 1.88 5.53 2.94T5 1.69 1.68 1.27 5.46 4.26T6 1.17 1.20 0.79 5.28 6.61T7 1.69 1.79 1.39 5.52 3.96T8 1.17 1.21 0.81 5.33 6.62T9 0.70 0.73 0.32 5.25 15.92R1 2.74e-04 2.98e-04 1.17e-04 1.49e-03 12.70R2 1.37e-04 1.42e-04 5.60e-05 7.OOe-04 12.63H 100.0 98.3 98.3 1.75 0.02F 50.0 70.8 27.9 386.24 13.73Hessian matrix calculated using the estimated parameter values. Either of the two setsof parameter estimates listed in Table 5.9 leads to confidence regions of similar shapeand size, though centered on different points in parameter space. Table 5.10 gives theaxes of the confidence region. Table 5.11 shows the relative contribution to the totalcoefficient of variation from each axis of the confidence region. Table 5.12 is the errorratio matrix for this problem.The longest axis of the confidence region is axis 13, which is 75 times longerthan any other axis. The valley in the multidimensional response surface is oriented100-I CD 9, -& 0 > >< CD Ci) 0 3 CD -4. CD -‘ C) 0 -4. a- CD D C.) CD CD 0 -4’0 3 C -4--o-‘ 3 CD -4. CD-o-‘ 0 0 CD 3-I 0Axis2.1e-08.2e-0C1.0e-02.2e-02.5e-02.7e-0:3.8e-05.6e-08.3e-02.Oe-012.2e-013.Oe-0122.30LengthsUnitU1U2U3U4U5U6U7U8U9U10U11U12U13VectorsT1-0.300.08-0.240.18-0.080.68-0.440.22-0.060.280.09-0.010.10T2-0.120.030.000.01-0.020.15-0.020.150.83-0.49-0.050.020.11T3-0.06-0.010.08-0.100.06-0.130.09-0.180.500.79-0.040.190.10T4-0.20-0.07-0.54-0.020.100.100.10-0.770.00-0.150.03-0040.10T5-0.13-0.14-0.09-0.32-0.070.390.770.26-0.120.03-0.060.1.00.14T6-0.11-0.240.17-0.700.450.05-0.370.00-0.09-0.09-0.080.060.22T7-0.15-0.06-0.660.040.25-0.48-0.020.480.000.050.050.050.13T8-0.11-0.35-0.08-0.24-0.81-0.21-0.21-0.02-0.04-0.03-0.020.190.22T9-0.10-0.660.240.540.230.020.04-0.03-0.06-0.050.000.280.30R10.200.500.020.020.020.00-0.02-0.05-0.13-0.130.110.640.50R20.070.120.000.11-0.04-0.04-0.020.00-0.060.05-0.76-0.380.49H0.86-0.29-0.31-0.060.010.24-0.090.030.120.060.000.010.00F0.07-0.010.13-0.04-0.04-0.060.070.000.030.040.63-0.530.53-1 CD 1 CD CD 0 0 -4.-‘ 0 -4.0 D -4.0 -4.CD C) -4’ 0 3 CD 0 Zr x Cl) 0 -4’ Zr CD C) 0 D -4’ a CD C) CD CD-D Cl) 0 a0ParametersU1U2U3U4U5U6U7U8U9U10U11U12U13T10.000.000.000.000.000.000.000.000.000.000.000.001.00T20.000.000.000.000.000.000.000.000.000.000.000.001.00T30.000.000.000.000.000.000.000.000.000.000.000.001.00T40.000.000.000.000.000.000.000.000.000.000.000.001.00T50.000.000.000.000.000.000.000.000.000.000.000.001.00T60.000.000.000.000.000.000.000.000.000.000.000.001.00T70.000.000.000.000.000.000.000.000.000.000.000.001.00T80.000.000.000.000.000.000.000.000.000.000.000.001.00T90.000.000.000.000.000.000.000.000.000.000.000.001.00R,0.000.000.000.000.000.000.000.000.000.000.000.001.00R20.000.000.000.000.000.000.000.000.000.000.000.001.00H0.010.020.030.000.000.110.030.010.370.380.000.020.01F0.000.000.000.000.000.000.000.000.000.000.000.001.00m -S -S 0 -S 1) 0 I0 (D 01 F’)0SpecifiedParametersEstimatedT1T2T3T4T5T6T7T8T9R1R2HFParametersT11.000.950.881.040.720.470.780.470.330.210.2240.20.22T21.051.000.921.090.760.490.820.500.330.220.2339.80.23T31.141.081.001.180.820.530.890.540.320.240.2539.80.25T40.960.920.851.000.700.450.750.460.320.200.2140.00.21T51.391.321.221.441.000.651.080.660.410.300.3049.20.30T62.132.031.872.211.541.001.661.010.590.460.4670.10.46T71.281.221.131.330.930.601.000.610.410.270.2849.20.28T82.112.011.862.191.520.991.651.000.590.450.4670.20.46T93.033.043.113.122.431.692.441.681.000.620.64101.0.57R14.614.394.064.793.332.163.602.191.631.001.00124.1.00R24.614.394.064.793.332.163.602.191.621.001.00124.1.00H0.000.000.000.000.000.000.000.000.000.000.001.00.00F4.614.394.064.793.332.163.602.191.741.001.00125.1.00in the direction of this axis, and is very long and flat, leading to high correlationsbetween any parameters that contribute to this axis. Examining the unit vector thatcorresponds to axis 13, all parameters except H have a component in the directionof axis 13. The axes of parameters F, R1 and R2 are most closely aligned with the axis13. Based on the contribution to the parameter CV matrix (Table 5.11), axis 13contributes nearly all the uncertainty to every parameter except H. If axis 13 could beshortened, the uncertainty in all parameters except H would be reduced.The linearized error ratio matrix (Table 5.12) shows that errors in parametersR1, R2, and F lead to small error ratios for the other parameters. Errors in thetransmissivity zones lead to larger error ratios for the other parameters. Errors in thehead boundary value lead to the largest error ratios.The confidence region for the parameter estimates can be interpreted asfollows. There is one axis that is much longer than the others, and this axiscontributes nearly all the uncertainty to all parameters except H. This axis is orientedmost closely with parameter axes F1, R1 and R2, and is perpendicular to parameteraxis H. This implies that parameter H is virtually independent of the other parameters,and can be estimated from the data even in this ill-conditioned problem.The unit vector associated with the longest parameter axis has some interestingproperties. The magnitudes of the components of the unit vector seem to fall into104three distinct groups. The first group contains the transmissivity parameters. Thecomponents of the unit vector range from 0.10 to 0.30, and seem to be inverselycorrelated with the magnitude of the estimated parameter. Since the confidence regionreflects the coefficient of variation, the smaller the value of the estimated parameter,the larger the coefficient of variation. The parameter estimates for T1 through T9 haveapproximately equal variances, so their coefficients of variation are inverselyproportional to their estimated values. In this confidence region, the longest axiscontributes nearly all the uncertainty to all parameters, so the components of the unitvector are closely correlated to the coefficients of variation. The second group ofcomponents contain the recharge and flux parameters. All three parameters havecomponents of the longest vector of similar magnitude, and thus similar coefficientsof variation. The variances of these parameters are quite different, but theircoefficients of variation are similar. This is the opposite of what was observed in thefirst set. The third group consists of the specified head parameter, which has nocontribution to the longest axis of the confidence region.5.6.2 Anticipated influence of prior information based on parameter spaceFrom the parameter space, three groups of parameters emerge. The headboundary value can be identified quite well using head data only. Prior information onparameter H will not reduce the ill-conditioning of the parameter set significantly.However, errors in the prior information for H would cause very large errors in the105other parameter estimates. Prior information for parameter H would be a poor choice,since it is not necessary, provides no advantage in identifying the other parameterestimates, and leads to large error ratios.The flux boundary value and the recharge parameters all have largecomponents of the unit vector associated with the long axis of the confidence region.These three parameters also have the smallest error ratios. Any of these threeparameters are good choices for prior information because they efficiently reduce theill-conditioning of the parameter set while reducing the influence of any errors in theprior information.The transmissivity parameters have smaller but still significant components ofthe unit vector associated with the longest axis of the confidence region. Priorinformation for any of these parameters will reduce the ill-conditioning of theparameter set significantly, but at the expense of potentially magnifying any errors inthe specified values of the parameters.5.6.3 Comparison of parameter space and modelled resultsBased on the analysis of the parameter space, the flux boundary and therecharge parameters should be the most efficient and responsible parameters for priorinformation. The transmissivity parameters are less efficient and responsible, while106prior information on the head boundary would be both inefficient and irresponsible.These results need to be checked against the actual modelled results.First, the issue of which parameters most efficiently stabilize the inverseproblem is examined. Two types of prior information are considered: prior informationwithout uncertainty and prior information with uncertainty. Table 5.13 shows the effectof specifying various parameters at their true values, with no uncertainty in the priorinformation. The component of the longest axis for each parameter is listed, along withthe average CV and the CN for the parameter set after specifying each parameter.Specifying those parameters with large components of the longest axis should resultin the most stable parameter sets with the smallest average coefficient of variation.Table 5.13 Calculated results from specifying some representative parameters attheir true valuesFull Specify Specify Specify Specify Specify SpecifyParameter T1 T5 T9 R1 F HSetComp. N/A 0.10 0.14 0.30 0.50 0.53 0.00of LmaxCN 9900.0 289.8 180.6 168.4 164.3 165.1 8944.0Average 5.510 0.181 0.113 0.102 0.097 0.099 5.490CvThe full parameter set (all 13 parameters) results in a very large CN and alarge average CV. Specifying parameter H does not significantly reduce either the CN107or average CV. Specifying any of the other parameters does significantly stabilize theproblem, because the dimension in parameter space in which all the parameters werecorrelated is removed. Table 5.13 shows that specifying parameter T9 reducesparameter uncertainty more than specifying parameters T1 or T5. Parameter T9 has alarger component of the longest axis than the other transmissivity parameters, sospecifying T9 should reduce parameter uncertainties more than any othertransmissivity parameters. The parameters with the largest components of the longestaxis are R1, R2, and F. Table 5.13 shows that specifying parameters R1 and F resultin the most stable parameter estimates with the lowest uncertainties. The behavior ofparameter R2 is similar to parameter R1. Specifying those parameters with largecomponents of the longest axis are generally the most efficient in stabilizing theparameter set.Table 5.14 illustrates the effect of prior information with uncertainty. Priorinformation is added to one parameter at a time, where the uncertainty in the priorinformation is ±25% of the prior value. The CN and average parameter CV arepresented for prior information on representative parameters. Prior information on R1and F resulted in the most stable parameter estimates with the lowest uncertaintiesof any parameters. If only transmissivity parameters are considered, prior informationon T9 resulted in the best parameter estimates. Prior information on parameters withthe largest elements of the unit vector associated with the longest axis of theconfidence region are most effective in reducing the average CV.108Table 5.14 Calculated results from prior information on various parametersAll Prior Prior Prior Prior Prior Prior HParameters T1 T5 T9 R1 FComp. of N/A 0.10 0.14 0.30 0.50 0.53 0.00LmaxCondition 9900.0 1181.9 973.6 413.2 274.0 295.0 8944.0NumberAverage 5.510 0.628 0.512 0.242 0.174 0.187 5.490CVTo determine for which parameters prior information will most responsiblystabilize the inverse problem, error ratios can be calculated. Prior information wasincluded one parameter at a time, where the prior value was 25% greater than thetrue value. Table 5.15 gives the calculated error ratios, calculated as the ratio of thepercentage error in the estimated parameter to the percentage error in the prior value.The first column of Table 5.15 are the calculated error ratios for the full set ofestimated parameters, where the starting estimate is equal to the true value of theparameter. The other columns are the calculated error ratios when one parameter ata time is specified at a value 25% higher than its true value. Because the data haveerrors, the calculated error ratios do not equal the linearized error ratios from Table5.12. However, the calculated error ratios do follow the trend of the linearized errorratios estimated from the parameter space. For instance, when T1 is specified, theerror ratios for T1 through T5 are near one, T6 through T8 are between one and two,T9 is above two, and the recharge and flux parameters are much higher. The109calculated error ratios for specified recharge and flux parameters are much lower thanthe error ratios for specified transmissivities. There is one pattern that the linearizederror ratios did not predict, that is the fact that the errors in the flux boundary valueas a consequence of errors in specified values for other parameters are consistentylarger than errors in recharge. The linearized error ratios predict that these errors willbe similar in magnitude. The cause of this discrepance may be the difference betweenthe true and estimated parameter values. The estimate of the value of the fluxparameter is much farther from its true value than the estimates of the rechargeparameters from their true values, even though the coefficient of variation for all theseparameters is similar in magnitude. This error due to data uncertainty may carry thewhole way through the error ratio analysis. No error ratios for H were calculated,because the parameters set was so ill-conditioned when using only prior informationon H that parameter estimates could not be reliably estimated.Because the parameter set is ill-conditioned without prior information, the finalestimates do not depend on the uncertainty assigned to the prior information. Evenif the prior information is assigned an uncertainty of 100%, the final estimates areidentical to those using prior information without uncertainty. Table 5.15 is unchangedfor any reasonable uncertainty in the prior information.If prior information is provided for more than one parameter, the situationbecomes more complicated. The relationship between the location of the prior110-‘—CDCD-‘a CD -.‘ 0 -‘ 0 C,) a CD 0 F’) 01 CD CD C’)-o CD C.) - CD a-1 0) 0 CD 9, - 01Errorratiocluetoerrorinspecified_parametersParametersTrueEstimateErrorSpecifySpecifySpecifySpecifySpecifySpecifyValuesRatioT1T5T8T9R1FT1150.00144.00-0.031.000.890.560.340.19-0.02T2150.00109.00-0.250.370.660.330.11-0.04-0.24T3150.0086.90-0.440.630.490.16-0.07-0.21-0.42T4150.00168.000.090.92-0.100.680.460.310.10T550.0041.80-0.180.932.420.570.280.10-0.17T615.0013.80-0.121.461.590.970.570.27-0.10T750.0054.400.091.231.260.840.550.370.10T815.0014.07-0.091.491.611.000.570.27-0.07T95.004.70-0.152.522.711.670.970.53-0.11R12.74e-042.58e-04-0.237.118.003.932.011.00-0.16R21.37e-041.24e-04-0.396.547.503.561.720.76-0.35H100.0098.30-0.07-0.05-0.07-0.07-0.07-0.07-0.07F50.0061.400.9111.7011.616.303.812.501.00estimates in parameter space, the relative uncertainties of the prior estimates, and thelocation of the valley in parameter space determine the final parameter estimates. Forinstance, when prior information on F, R1 and R2 is provided, and all three priorestimates are 25% greater than their true value, the final parameter estimate for R1and R2 are approximately 15% greater than their true values, and the final estimatefor F is approximately 35% greater than its true value. In another example, when theprior estimate of F is 25% less than its true value, and the prior estimates of R1 andR2 are 25% greater than their true values, then all three final estimates are within 5%of the true parameter values. The errors in the prior estimates are averaged, and thefinal estimates do not reflect the errors in the prior estimates.5.6.4 Discussion of results for multi-parameter systemIn order to responsibly stabilize this problem, the parameters F, R1 and R2 werefound to be the most efficient and responsible parameters to select for priorinformation. Practically, these three parameters may be the most difficult to measurein order to obtain the prior information. It is difficult to estimate the value of a fluxboundary independent of the model. Recharge is also a difficult parameter to estimatewithout the use of the model unless a detailed field measurement program isimplemented. However, errors in the prior estimates of these parameters translate intorelatively small errors in the estimates of transmissivity parameters. Prior estimateson these three param,even if they are in error or have large uncertainties, may112be better for parameter estimation than prior information in the transmissivityparameters. A systematic underestimation of prior estimates on transmissivityparameters may lead to large errors in the estimation of the recharge and fluxparameters.Should prior information on parameters be used if it is highly uncertain or hasthe potential for significant errors? Three factors need to be considered: theconditioning of the parameter set without prior information, the error ratio of theparameters with prior information, and the distribution of the prior information. If theparameter set is ill-conditioned without the prior information, then the prior informationwill constrain the model parameter estimates. If the error ratio is large for theparameters with prior information, then errors in the prior values may lead to largeerrors in the estimated values of the other parameters. If the prior information has thepotential to systematically over or underestimate the true parameter values, then theprior information may lead to errors in parameter estimates. If the prior information israndomly distributed about the true value of the parameter estimate, then the priorinformation might not lead to significant errors in the final parameter estimates.However, the assumption that prior information is randomly distributed about the trueparameter value is rarely fulfilled in practice.In this synthetic system, the full parameter set is ill-conditioned without priorinformation. Parameters R1, R2 and F lead to small error ratios, but these are difficult113parameters to measure independent of the model. The transmissivity parameters leadto larger error ratios, but are easier to measure in the field. However, it is common forprior information on transmissivity to be significantly biased, because of scalingdifferences between the measurement of transmissivity in the field and the modeltransmissivity. One solution to this dilemma is to provide prior information onparameter T9. This parameter leads to the smallest error ratios of any of thetransmissivity parameters. Errors in this parameter value will not be magnified asmuch as errors in any of the other transmissivity parameters.If at this point the parameter set was still unstable, or a further reduction inparameter uncertainty was required, the same analysis could be carried out to selecta second parameter for prior information. The reduced parameter space (12parameters now) would be analyzed, and a parameter selected. This process is easyto accomplish, since the parameter estimation in the 12 parameter system has alreadybeen done. The scaled Hessian matrix at the parameter estimates would be used toconstruct the three tables given above for the 12 parameter system. This parameterspace would be analyzed, and a second parameter selected.5.7 Summary and conclusionsParameter estimation problems are often non-identifiable or unstable using onlyhead data. Prior information on some of the parameters may be used to stabilize the114parameter set for the purpose of parameter estimation. Guidelines have beendeveloped for efficient and responsible use of prior information in parameterestimation.An analysis of the parameter space is used to select parameters for priorinformation. The most efficient parameters for prior information are those with thelargest element of the unit vector (eigenvector) associated with the longest axis of theparameter confidence region. The axes of these parameters in parameter space aremost closely aligned with the longest axis of the confidence region. To minimize theeffects of error in prior estimates, the most responsible parameters for priorinformation are those which lead to the smallest error ratios. The linearized error ratioscan be calculated from the axes of the parameter confidence region. The parametersthat lead to the smallest error ratios are often the same parameters that have thelargest element of the unit vector associated with the longest axis of the parameterconfidence region.The parameter space analysis was applied to a multi-parameter system withgood results. The parameter space was used to gain an understanding of theinteraction of the model parameters and the data set. The parameters for which priorinformation most efficiently and responsibly stabilized the parameter set were correctlyidentified using the parameter space. The linearized error ratios followed the patternobserved in the true error ratios. The parameter space analysis may lead to selecting115parameters about which reliable prior information is difficult to obtain. The modellermust make a judgement regarding the reliability of the prior information against thepotential effects of errors in the prior estimate. The parameter space analysis givesthe modeller a tool to judge the potential effects of errors in the prior estimates.116ECDCDC\J‘Ici)ci)C)ci)(1)Transmissivity Zone 1T, = 20 m2ldayTransmissivity Zone 2T2 = 200 m2IdLI Boundary of Recharge ZoneR = 0.0004 mId13 Sample LocationFigure 5.1 Five parameter synthetic flow system.1172CDCDCDCDSpecified Head = 100 m0.85 0.88 0.91 0.94 0.97 1 .00 1 .0.3 1 .06 1 .09 1 .1 2 1 .1 51.15 I 11111111111111 I I ii III it 1.15E1.00 / ///// / E 1.000.97 - - 0.97085 0.850.85 0.88 0.91 0.94 0.97 1.00 1.0.3 1.06 1.09 1.12 1.15T1Figure 5.2 Response surface for T1-R parameter set using hydraulic head data1180.85 0.88 0.91 0.94 0.97 1.00. 1.03 1.06 1.09 1.12 1.151.15 1.151 .06R 1.001.120.85 0.850.85 0.88 0.91 0.94 0.97 1.00 1.03 1.06 1.09 1.12 1.15T1Figure 5.3 Response surtace for Ti -R parameter set using prior information on R.1.121 .091.03liii III II18.4 18.4 =9.2 9.2—4.6 4.6-2.3 2.3—1:1 1- 2.3 2.34.6 4.6- 9.2 9.2- 18.4 18.4I I I I I I I I I I I I I I I I I I I I I0.970.940.910.881.091.061.031.000.970.940.910.881190.85 0.88 0.91 0.94 0.97 1.00 1.03 1.06 1.09 1.12 1.151.15—......j———— I I I I I I I I I I I I I I I I I I I I 1 1.151.12- 1.121.09—1.091.06 / I - 1.06/ •/i -/ i/I II4- —1.03 /- 1.031.:::: 7 L\0/ ::0.85 1It 0.850.85 0.88 0.91 0.94 0.97 1.00 1.03 1.06 1.09 1.12 1.15T1Figure 5.4 Response surface for T1-R parameter set using head and priorinformation on R1200.95 0.96 0.97 0.98 0.99 1.00 1.01 1.02 1.03 1.04 1.05i.O4 - 1 .041.05_I II Ijl I,,/I II 11111051.03 - 1 .031 .02- (01.01 - 1.011 02 . .T2 1.00 - 1 .000.99 - 0.990.98 - 0.980.97 - 0.97CD CJ0.950.96 0.960.950.95 0.96 0.97 0.98 0.99 1.00 1.01 1.02 1.03 1.04 1.05T1Figure 5.5 Response. surface for T1-T2 parameter set using hydraulic head data1211.50.5• +2.5% Parameter estimates witherrors in head B.C. parameter• +25% Parameter estimates witherrors in recharge parameterFigure 5.6 Effect of errors in head and recharge parameter values on estimates ofT1 and T2T21.00.51.0T11.5122T2T,Figure 5.7 Orientation of longest axis ofT1-T2H confidence ellipsoid123T2Figure 5.8 Orientation of longest axis of T1-T2R confidence ellispoid124A(Y-Y0)Figure 5.9 Generic confidence ellipse to show calculation of error ratio for parameterb2 due to an error in parameter b1(X,Y)b1A’125AError in b2 due to(X-x,,) error in b1(a) Poorly-conditioned confidence ellipse.b2Lon axis ofcon ide::e elI7./Ic(X-X0)AError in b2 due to(X-X0) error in b1A’(b) Well-conditioned confidence ellipse.Figure 5.10 Influence of shape of confidence ellipse on error ratio for (a) well-conditioned and (b) poorly-conditioned confidence regionsofb1260.95 0.96 0.97 0.98 0.99 1.00 1.01 1.02 1.03 1.04 1.051.05 I I I I I I I I :1051.04-- 1.04N) Q Ci--1.03- 1.031.02- 1.02- O CD(N1.01 - 1.01T2 1.00 - 1.000.99 - o r 0.99N) 0)0.98- 0.980.97 - 0.97- CD0.96- 0.960.95 II I I I I I Q950.95 0.96 0.97 0.98 0.99 1.00 1.01 1.02 1.03 1.04 1.05T1Figure 5.11 Response surface for T1-T2 parameter set using incorrect priorinformation on T21270.85 0.88 0.91 0.94 0.97 1.00 1.03 1.06 1.09 1.12 1.151.15 I I I I I I /i ,i,,/’ii ) I 1 I 1.15• ol?4(39 / //// )/ 1.091.06 / /// J/ 1.061.03 / /// /// 1.031.00 / //L/ / 1.000.97 . t-1/ / 0.970.94 0.940.91 / 0.9100.88 0.880.85 I I I I I I L I I I I /f” I I I I I I I I I I I 0850.85 0.88 0.91 0.94 0.97 1.00 1.03 1.06 1.09 1.12 1.15T1Figure 5.12 Response surface for T1-R parameter set using head and incorrect priorinformation on T11280.95 0.96 0.97 0.98 0.99 1.00 1.01 1.02 1.03 1.04 1.051.05 I III I/I I,,/ 111111 I I I L105104 0- 1.041.03- 1.03-- 1.021.02I /-uc311.01 1.01T2 1.00 ( 0.991.00j()10.99 0 cc)0.98 0.980.97 0.97\I I I II\ 0.960.96 I I I0.95 0.950.95 0.96 0.97 0.98 0.99 1.00 1.01 1.02 1.03 1.04 1.05T1Figure 5.13 Response surface for T1-T2 parameter set using head and incorrect priorinformation on parameter T2129EU)c4.JII><U--DC)C.)C)Cl)600CmFigure 5.14 Multi-parameter synthetic flow systemtransmissivity and recharge zones130showing arrangement ofNT7 T8 T950 15 5T4 T5 T6150 50 15T1 T2 T3150 150 150Specified Head = 100 mR1 = 0.000274 m/dR2 = 0.000137 mIdTransmissivity Zones(m2/d)Recharge ZonesCHAPTER 6. PARAMETER SPACE METHODS IN JOINT PARAMETERESTI MATIONParameter estimation for hydrogeological models using hydraulic head data isfrequently plagued by difficulties due to the illposed nature of the inverse problem.In Chapter 5, guidelines were developed for using prior information to overcome thesedifficulties in an efficient and responsible manner. Joint parameter estimation providesanother method for overcoming these difficulties. Incorporating a second set ofobservation data of a different data type through joint parameter estimation may assistin stabilizing the inverse problem.Joint parameter estimation extends the concept of single state parameterestimation by incorporating more than one set of observation data in the parameterestimation process. In this thesis, the two data sets are hydraulic head data and tracerconcentration data. The tracer concentration data is used to help estimate the flowparameters. Transport parameters are not estimated, and need to be specified.However, effects of errors in the specified transport parameters can be evaluated. Thegoverning equations and boundary conditions for the simulation of flow and masstransport have been introduced in Chapter 1. This chapter examines joint parameterestimation in the context of parameter space analysis.Three issues in joint parameter estimation are examined in this chapter. First,131parameter space analysis is used to show how the two data sets interact to yieldimproved parameter estimates. Second, parameter space information can be exploitedto determine a priori if a future data set will improve parameter estimates significantly,before the second data set is collected. The third issue explored is the weighting ofthe different data sets in joint estimation. Four methods of weighting data sets for thepurposes of minimizing parameter uncertainty with the available data are presentedand evaluated.6.1 Parameter space analysis of multiple data setsParameter space analysis is useful for understanding how two data sets interactto produce an improved set of parameter estimates. This improvement is manifestedin a better conditioned, more stable set of parameter estimates with a smallerparameter confidence region. The hydraulic head data set defines a response surface.If, the model parameters are poorly conditioned, the response surface will contain avalley with a flat bottom. The response surface defines the parameter confidenceregion for the head data. A second data set, the tracer concentration data, defines anew response surface. The total response surface used in joint estimation is the sumof the individual response surfaces. If the tracer concentration response surface hasa shape and orientation that is significantly different from the head response surface,then the total response surface will be better conditioned than the initial headresponse surface. If the tracer concentration response surface has a shape and132orientation that is similar to the head response surface, then the parameter confidenceregion will not be reduced significantly.The confidence regions defined by the response surfaces can be approximatedby the linearized confidence ellipsoids (Equation 4.4). The difference in orientation ofthe confidence regions for the two data sets can be expressed in terms of theirconfidence ellipsoids. The angle between the longest axes of the two confidenceellipsoids can be used to approximate the difference between the two confidenceregions. This angle, called the angle of interaction (j3), can be calculated using:cos (f3) = Uhi. U1where Uh is the unit vector corresponding to the longest axis of the confidenceellipsoid using the head data set, and U is the unit vector corresponding to thelongest axis of the confidence ellipsoid using concentration data set. A small angle ofinteraction represents confidence regions with similar orientations, while a large angleof interaction represents confidence regions with very different orientations. The largerthe angle of interaction, the greater the potential increase in stability and decrease inthe uncertainty of the parameter estimates using the joint data set over the initial dataset.6.1.1 Demonstration using a two parameter setThe five parameter synthetic flow system introduced in Chapter 5 is used to133demonstrate response surface analysis for multiple data sets. Tracer transport issimulated within this flow system. A specified concentration is introduced along theinflow head boundary, and the tracer concentrations are measured at the samelocations as the head data after 500 days. For this simulation, all units are assignedan effective porosity of 0.2, and the specified concentration parameter values areczL=100 m, T = 10 m, and c1 = 1.0. Figure 6.1 shoWs the concentration distributionthroughout the flow system at 500 days. The leading edge of the plume has movedmost of the way through the flow system. Some of the concentration samples, in thecenter of the lower permeability zone, are not sensitive to changes in the parametervalues. The simulated measurements of tracer concentration are corrupted withGaussian distributed errors with an uncertainty of 5.0%.Response surtaces are used to show how information from the two data setsinteract to produce smaller confidence regions in joint parameter estimation. Only twoparameters are estimated in order to visualize the response surfaces. Responsesurfaces and confidence regions are discussed for two sets of parameters: (1) Set A:transmissivity zone 1 (T1) and recharge zone (R); and (2) Set B: transmissivity zone1 (T1) and transmissivity zone 2 (T2). These two parameter sets provide goodexamples of the influence of a second data set on the uncertainty and stability of theparameter estimates.Table 6.1 lists the parameter estimates and uncertainties, along with the134Table 6.1 Parameter estimates and uncertainties for parameter sets A and BHead Data Concentration CombinedData DataT1 Estimate 1.320 1.326 1.320R Estimate 0.00042 0.00042 0.00042Std (T1) 0.0498 0.0956 0.0430Std (R) 4.4e-05 9.4e-05 3.9é-05CV(T1) 0.0377 0.0619 - 0.0345CV (R) 0.1064 0.21 20 0.0988Average CV 0.0721 0.1370 0.0667T1 Estimate 1.303 1.303 1.303T2 Estimate 2.284 2.297 2.297Std (T1) 0.0102 0.0143 0.0079Std (T2) 0.0468 0.0089 0.0078CV (T1) 0.0079 0.0110 0.0061CV (T2) 0.0204 0.0043 0.0034Average CV 0.01 41 0.0076 0.0047coefficient of variation, for the two parameter sets. The parameter estimates for T aregiven in log-transformed values, and the parameter estimates for R in units of m/d.Equal weights are assigned to each data set for joint estimation. Table 6.1 containsthe estimates for both the individual data sets and the joint data sets.An averagecoefficient of variation is included to allow easy comparison between the data sets.This table shows that for parameter set A, the additional concentration data does notreduce the parameter uncertainties significantly, as measured by the coefficient ofvariation. The average CV for the combined data set is only slightly better than for the135head data set. For parameter set B, the additional concentration data does reduce theparameter uncertainties significantly. The average CV for the combined data set ismuch smaller than for either the head or concentration data sets. An analysis of theresponse surfaces shows why the two parameter sets behave differently under jointparameter estimation. -The response surfaces for parameter sets A and B are shown in Figure 6.2 and6.3. Each figure has three parts: (a) is the response surface using only head data, (b)is the response surface using only concentration data, and (c) is the total responsesurface. The 68% joint confidence region is defined by the area inside the 2.3 contourof the response surface.The response surface for parameter set A (Figure 6.2a) using only head datacontains a long, narrow valley with a nearly flat bottom. The estimates of T1 and R arehighly correlated, and the parameter set is ill-conditioned. The confidence region usingconcentration data only (Figure 6.2b) is similar in shape and orientation to that usinghead data, but is larger in size. The confidence region for the combined data sets(Figure 6.2c) is slightly smaller than that for head or concentration data alone, and isvery similar in shape and orientation. The angle of interaction calculated from thelinearized confidence ellipsoids is 3 degrees, so the orientation of the confidenceellipsoids are very similar. Since the two confidence regions are very similar, additionof the second data set does not significantly increase the stability nor reduce the136uncertainty of the parameter estimates.The response surface using head data only for parameter set B (Figure 6.3a)is better conditioned than that for parameter set A. A definite minimum Is observed,bat it. stIll contains a valley oriented almost parallel to the T2 axis. The confidenceregion using concentration data (Figure 6.3b) Is oriented nearly parallel to the T1 ads,and is smaller than the confidence region using head data. The angle of Interactionfor this parameter set Is 77 degrees, confirming that the orientations of the twoconfidence regions are significantly dIfferent The confidence regIon for the combineddata set (Figure 6.3c) Is much smaller than the confidence regions for the lndMdualdata sets. The parameter estimates are stable, and the uncertainty In the parameterestimates has decreased substantially The data sets can be combIned to yield stable,well-defined parameter estimates.For parameter sets A and B, the addition of concentration data led to verydifferent behaviors. With parameter set A, the concei,tiaflon data did not reduce theuncertainty of the parameter estImates sIgnificantly. With parameter set B, thecombined data sets resulted In stable parameter estimates with low uncertainties. Theangle of interactIon reflects the dIfference in orIentatIon of the confidence regions.it Is often difficult to estimate both transmlsslvity and recharge together [eg.Woodbury et aL, 1987]. This analysIs shows that when estimating T1 and R137simultaneously, the head data and the concentration data contain similar informationabout the model parameters. Both data sets are very sensitive to the value ofparameter T1, and not very sensitive to the value of parameter R. The joint data setdoes a poor job in estimating the two parameters. Another example, estimating T2 andR simultaneously, is not presented here in detail but shows a different result. For theparameter set T2 and R, the head data and concentration data contain differentinformation about the parameter values. The response surfaces are different in shapeand orientation, with an angle of interaction of 39°, and the joint parameter estimatesare more stable than either of the single state estimates. There seems to be no clearphysical explanation as to why, in one case, the head and concentration data containsimilar information about the parameters, and in the second case the head andconcentration data contain different information about the parameters.6.1.2 Multiple parameter dimensionsIn most parameter estimation problems, more than two parameters areestimated. In these cases, the response surfaces cannot be visualized. The responsesurface in the form of the parameter confidence ellipsoid can be approximated. Thelongest axes of the parameter confidence ellipsoid correspond to the directions ofmaximum uncertainty in parameter space.As with the two parameter problem, the confidence ellipsoids can be calculated138for both data sets. The topologies of the ellipsoids for the two data sets can becompared by comparing the lengths and directions of the longest axes. If the longestaxis of the confidence regions for the two data sets are similarly oriented and similarin size, the second data set will not reduce the parameter uncertainties or instabilitiessignificantly. If the longest axes are similarly oriented but of different sizes, the jointparameter estimates will be similar to those of the data set with the smaller confidenceellipsoid. If the longest axes are oriented in significantly different directions, the seconddata set has the potential to significantly reduce the parameter uncertainties. However,simply using the orientation as measured by the angle of interaction between thelongest axes can be misleading in multiparameter estimation problems. Several axesof the confidence ellipsoid may contribute to the ill-conditioning or large uncertaintiesin the parameter estimates, and these other axes are not taken into account using theangle of interaction.It may be helpful to identify the largest axes for each data set by examining theeigenspace, and then calculate the angle of interaction between each of these axes.For instance, if each confidence ellipsoid had two axes that were much longer thanthe others, a total of four angles of interaction would be calculated. If all angles ofinteraction were large, then the ellipsoids would be significantly different. If one or twoof the angles were small, then the combined ellipsoid would still contain a directionwhere the parameter confidence had not been reduced very much. The parameterestimates may not improve significantly.139A second method of comparing the confidence ellipsoids is to compare theeigenspaces for the two data sets. The magnitude of all axes of the confidenceellipsoid for the head data set should be examined to see how many contributesignificantly to the uncertainty of the parameter estimates. The relative contribution tothe uncertainty of each parameter from each of these long axes can be calculated.The directions of the largest axes for the head data set are then compared to the longaxes for the concentration data set. This method has the advantage of a morethorough analysis of the parameter space, leading to a better understanding of themodel. However, it has the disadvantage of not producing a single number that allowsthe comparison between the two confidence regions. These methods are used insection 6.4, when a multi-parameter example is examined.6.1.3 Multiple sampling periodsThe above examples used tracer concentration data sampled at a single pointin time for parameter estimation. These samples provide a single snapshot in time ofthe concentration distribution. Numerical experiments were also conducted for thetracer concentration data sampled at multiple time intervals, at 100, 500, 1000 and1500 days, as the tracer moved through the system.The response surfaces with the concentration data sampled at multiple timeintervals are similar in shape and orientation to the response surfaces sampled at a140single snapshot in time. There are two differences between the response surfacespresented for a single sampling period and the response surfaces for multiplesampling periods. First, the response surfaces for multiple sampling periods result insmaller confidence regions, because more data exists. Second, the confidence regionsfor the data from the multiple sampling intervals are generally better conditioned thanthe confidence region for the single sampling period. For the trial with a singlesnapshot in time, some of the sampling locations are insensitive to changes in theparameter values. For the trial with multiple sampling periods, all of the samplinglocations are sensitive to changes in parameter values for at least one samplingperiod. Since the response surface for each sample location is slightly different, thetotal response surfaces for the concentration data are somewhat better conditionedwhen multiple sampling periods are used. Numerical experiments were also conductedfor a decaying environmental tracer, sampled at its steady state concentrationdistribution. The response surfaces for the environmental tracer are very similar to theresponse surfaces for a single snapshot in time.6.2 Predicting the usefulness of a second data setThe use of response surfaces described above can be extended to predict howadditional data will influence the uncertainty of the parameter estimates before thedata is collected. The basis for this prediction is an analysis of the parameter spacefor both the actual data set and a proposed future data set. The initial data set will141produce a set of parameter estimates along with the parameter confidence region inthe form of a response surface. The general shape of the response surface for thefuture data can be calculated even without knowing the actual values of the futuredata, although the position of the confidence region in parameter space cannot bedetermined until the data is collected. The superposition of the two response surfaceswill produce the total response surface for both the initial data and the proposed futuredata, from which the size and shape of the confidence region can be calculated. If theresponse surfaces for the initial and future data sets are similar, the future data setwill not significantly reduce the uncertainty in the parameter estimates. If the responsesurfaces are quite different, the addition of the future data set will improve theparameter estimates significantly.Even though the actual values of the future data are unknown, and the valueof the resulting parameter estimates are unknown, the general shape and size of theparameter confidence region can be calculated. The topology of the response surfaces(and thus the confidence regions) depend on how the data interact with the modelparameters. The location of the minimum in the response surface depends on theactual data values, but the overall shape of the response surface depends on theinteraction of the data with the model parameters. If the true data values are unknown,then the location of the minimum in the response surface is unknown. However, aslong as the physical location of the data samples are known, the general topology ofthe response surface can be calculated. The size and shape of the confidence region142for the parameter estimates can then be calculated.The general procedure is:1. Design a model and estimate parameters based on the initial data set.2. Use the model and parameter estimates to simulate a proposed future data set.3. Calculate the confidence regi6n for the future data set.4. Compare the confidence region for the initial data set to the confidence regionfor the future data set, using either the angle of interaction or analysis ofconfidence ellipsoids.5. Evaluate the potential reduction in parameter uncertainty based on thedifference between confidence regions for the two data sets.6.2.1 Demonstration using the T-T2 parameter setTo illustrate the above method, the parameter setT1-T2can be used. The initialdata set consists of hydraulic head data at the 15 observation locations, with a 1meter standard error. The proposed future data set includes 15 concentration datataken at 500 days after the source is introduced, It is assumed that the futureconcentration data have a 5% standard error due to measurement uncertainty. Thisexample is similar to the one presented above, except for the assumption that theconcentration data has not been collected before the process of parameter estimation.143Estimates for T1 and T2 are obtained using the initial hydraulic head data set(from Table 6.1). The response surface for the hydraulic head data has beenpresented in Figure 6.3a. These estimates of T1 and T2 are used in a forward modelto calculate the synthetic concentration distribution, and the concentrations aresampled at the observation locations. These synthetic concentration data, with noadded error, are used to estimate the parameters again. This time, the parameterestimates are not the major focus, since the estimates are the same as the estimatesusing the initial head data set. The parameter values obtained from the head data setwere used to simulate the future concentration data, so these parameter values arerecovered when using the synthetic concentration data to estimate the parameters.Instead, the parameter confidence region is the main concern. The response surfaceis calculated for the synthetic concentration data (Figure 6.4a). Comparison of theresponse surface for the synthetic concentration data to the response surface usinghead data (Figure 6.3a) shows that the orientation of the valley in the responsesurfaces are quite different. The joint confidence region (Figure 6.4b) is much smallerthan either region considered separately. The point estimates are the same as for thehead data alone.The synthetic future data set will not have the same values as the actual futuredata set. Therefore, the parameter estimates based on the synthetic data set will notbe the same as for the actual future data set. The general shape of the responsesurface for the actual future data set and the synthetic future data set will be similar,144although the location of the minimum are different. For example, the actual future dataset may have a response surface shown in Figure 6.4c. (This response surface wasgenerated by adding Gaussian noise to the synthetic concentration data, and thusrepresents a possible concentration data set.) Note that the shape of the responsesurface is similar to the synthetic future data set (Figure 6.4a), but the location of theminimum is different. The combined response surface for the head data set and theactual future data set is shown in Figure 6.4d. It is similar in shape to the combinedresponse surface for the initial data set and the synthetic future data set (Figure 6.4b),but the location of the minimum is different.In the above example, the future data set would significantly improve theparameter estimates. However, for an example using the T1-R parameter set, thefuture data would not improve the parameter estimates significantly. The responsesurfaces for the initial and future data sets are quite similar, so the future data will notsignificantly reduce the parameter ill-conditioning or uncertainty. It may not be worththe expense of collecting additional data in this case.6.2.2 Discussion of worth of future dataThis section has presented a method for evaluating the worth of future data forthe purpose of reducing parameter uncertainty. The size of the confidence region forthe combined present data set and the proposed future data set can be evaluated.145These confidence regions are compared to determine if the future data willsignificantly impact the parameter uncertainties. The actual parameter estimatescannot be calculated, since the true values of the future data are unknown. Differentsets of future data can be evaluated with respect to reducing parameter uncertainty.There are two major assumptions implicit in this method. First, the model forthe initial data set and the model for the future data set are the same. The modelparameters and their uncertainties are being compared, so the model, including thezonation and boundary conditions must remain unchanged. The initial data set shouldallow the construction of a model and estimation of reasonable parameter estimates.The future data are being evaluated on their ability to improve the parameterestimates.The second assumption is that the method is based on a linear analysis ofparameter uncertainty. If the model is highly non-linear with respect to the parameters,then the confidence regions become curved and non-symmetric. In these cases, theshape of the region depends on the value of the parameter estimates. If the model islinear with respect to the parameters, then the shape of the ellipsoid is independentof the value of the parameter estimates. Evaluating the worth of future data usingparameter space analysis relies on the fact that the response surface does notchange significantly when the parameter values change. Even in the case where themodel is highly-nonlinear in its parameters, this method may provide a first146approximation to the worth of the future data in reducing parameter uncertainty.The above analysis allows the worth of the future data to be quantified basedon its ability to reduce parameter uncertainties. The reduced parameter uncertaintieswill result in a model that has more accurate predictive capabilities. The benefits ofthe data must be weighed against the cost of collecting the data. A decision analysisframework may be used to determine the worth of the data based on the ultimatepurpose of the model [Freeze eta!., 1990]. The future data is only worthwhile if thebenefits produced by the future data exceeds the cost of the data collection.6.3 Weighting data sets in joint parameter estimationThe key difference between joint parameter estimation and single stateparameter estimation is the presence of multiple data sets. With multiple data sets,a decision must be made regarding the importance of each of the individual data sets.Each data set contains different, possibly unique, information on model parametervalues and uncertainties, yet some data sets contain more information than others.The importance of each data set is manifested as a set of weights in the parameterestimation process, as in equation (3.4). The data set containing more, or better,information should be weighted more heavily than data set containing less information.However, even data sets with poor quality information may be important since theymay contain unique information not found in the other data sets. It is important to147assign an appropriate set of weights to each data set in order to obtain the best setof parameter estimates possible. A well designed set of weights will extract themaximum information from each data set. A poorly designed set of weights mightignore essential information contained in some of the data sets. The use of parameterspace information allows the modeller to design a set of weights during the parameterestimation process.In this section, several methods designed to select an appropriate set ofweights in a joint parameter estimation problem are investigated. An analysis of dataresiduals is the traditional method for weighting of the data sets in joint parameterestimation [Galley eta!., 1992]. Using this method, a set of weights is selected sothat the residuals for each data set have approximately equal variance. In addition,three parameter space methods are proposed to select a set of weights. Theseparameter space methods select a set of weights based on an analysis of theparameter space of the individual data sets. Weights can be selected to satisfy anyof the following goals:1. Maximize parameter stability for the final parameter estimates.2. Minimize the total uncertainty of the final parameter estimates.3. Minimize the longest axes of the parameter confidence region.Each of the parameter space approaches are designed to accomplish different goals,and may lead to different sets of weights. All of the approaches are based on ananalysis of the parameter space. The final set of weights is selected to extract the148most information from each data set with respect to the goal that is required.6.3.1 Analysis of data residualsThe analysis of data residuals method uses the data residuals to determine, therelative weights of the data sets. This method was first proposed for using bothhydraulic head and prior information on parameters [Neuman and Yakowitz, 1979]. Ithas also been used with joint head-thermal data [Woodbury and Smith, 1988], andjoint head-concentration observations [Sun and Yeh, 1990; Galley eta!., 1992]. Theweights are assigned so that the variance of the weighted data residuals for each datatype are approximately equal.The following is an outline of the above approach for two data sets. For thejoint parameter estimation problem using head and concentration data the objectivefunction is:5(b) Wh(Yh—f(b)) V,,1 (Y,,—f(b)) + W(Y_f(b))T V (Y—f(b)) (6.1)where the h subscript denotes those terms relating to head data and the c subscriptdenotes those terms relating to concentration data. The Vh and V matrices containthe relative reliabilities of the data within each data set. The Wh and w terms areweights for each portion of the objective function, or each data set. The weights Whand w are to be determined. The general procedure for determining the weights is asfollows (although it varies slightly from author to author):1491. Minimize the objective function S(b) with Wh = w = 1.02. At minimum S(b), calculate variance of weighted data residuals for each dataset separately.WhSh W SSh2 = s2 = c cwhere Sh and S, are the minimum of the objective function using each data setindividually, and nh and n are the number of data in each data set.3. Compute the ratio of the variances.s2 hRatio = —4. Adjust the weights so that the ratio of the variances is equal to one, andestimate the parameters with the new set of weights.5. Repeat steps 2,3,and 4 until the ratio of the variances stabilizes [Neuman andYakowitz, 1979] or is approximately equal to one [Galley et a!., 1992]. Theiterative steps are necessary because the data residuals may change withdifferent sets of weights.6. Calculate the final parameter estimates and uncertainties.The key point for the analysis of data residuals method is that the weights are chosenso that the variances of the weighted residuals for each data set are approximatelyequal. In the examples which follow, the criterion which chooses weights based on theanalysis of data residuals will be termed the RESID criterion.1506.3.2 Methods based on analysis of parameter space.The parameter space contains information about the parameters beingestimated. This information can be used to select the weights in joint parameterestimation problems. Three methodshave been developed to choose the weights foreach data set. Each method emphasizes a different aspect for the final parameterestimates. The three goals for the final parameter estimates are:1. Parameter estimates with maximum parameter stability.2. Parameter estimates with minimum total uncertainty.3. Parameter estimates with the largest uncertainties minimized.The sections 6.3.2.1 through 6.3.2.3 outline the concepts for each of the weightingmethods, and section 6.3.2.4 outlines the procedural details for calculating the weightsfor all parameter space methods.6.3.2.1 Maximum parameter stabilityStable parameter estimates have the characteristic that uncertainties in the datavalues do not influence the value of the parameter estimates greatly. When parameterestimates are unstable, small errors in the data values can cause large errors in thevalues of the parameter estimates.Each data set has a unique response surface. Often these response surfaces151for individual data sets are somewhat unstable, containing a narrow valley. When theorientation of the valleys are different for each data set (measured using angle ofinteraction), the data sets can be combined to yield more stable parameter estimates.A specific combination of the data sets will produce parameter estimates with themaximum parameter stability.In order to choose weights to obtain parameter estimates with maximumstability, the condition number (CN) of the scaled Hessian matrix is minimized. TheCN for each data set is the ratio of the longest to the shortest axis of the parameterconfidence ellipsoid. Using joint parameter estimation, the shape of the confidenceellipsoid will change based on the weighting of each data set in the joint problem. Theweights that minimize the CN of the joint problem are to be selected.The advantage of this weighting scheme is that the parameter estimates areas stable as possible, given the model being calibrated and the data available. Itsmajor disadvantage is that the method may assign greater weights to the data set witha larger confidence region just to produce circular confidence regions. The methodmay produce stable parameter estimates with larger uncertainties than other weightingschemes.1526.3.2.2. Minimum total parameter uncertaintyThe overall parameter uncertainty can be measured by calculating the volumeof the parameter confidence region. To minimize the overall parameter uncertainty,the volume of this region is minimized [Sun and Yeh, 1985]. The volume of the regionis calculated using the axes of the parameter confidence ellipsoid, with the volumeproportional to the product of the lengths of the axes. During joint parameterestimation, the shape and size of the confidence ellipsoid will change based on theweighting of each data set in the joint problem. The weights that minimize the volumeof the confidence region for the joint problem are identified.The advantage of this method is that the final set of parameter estimates havea confidence region with the smallest possible volume. Having confidence region witha small volume is good, and this type of criterion is often used in experimental designfor the purpose of discriminating among competing models [Hill, 1978]. However, thesmallest volume confidence region does not necessarily lead to individual parameterestimates with minimum uncertainty. As an example, consider two confidence regionswith the same volume, Figure 6.5a and 6.5b. It is apparent that the confidence regionin Figure 6.5a leads to individual parameter estimates with much larger uncertaintiesthan those in Figure 6.5b. The parameter estimates represented by Figure 6.5a aremore unstable than those represented by Figure 6.5b.153Selecting weights which minimize the volume of the parameter confidenceregion often tends to reduce the smallest axes of the parameter confidence ellipsoidbefore reducing the largest axes. Reducing the smallest axes of the confidenceellipsoid reduces the volume much faster than reducing the largest axes. The resultof minimizing the volume of the parameter confidence ellipsoid is that the stability ofthe estimated parameters is often decreased.6.3.2.3 Minimizing longest axes of parameter confidence region.The above weighting methods minimize either the volume of the confidenceregion or maximize the parameter stability, often at the expense of the other criterion.A third weighting method, selecting weights that minimize the length of the longestaxes of the parameter confidence ellipsoid, strikes a compromise. When the longestaxes of the parameter confidence region are minimized, the parameter estimates tendto become more stable. The individual parameter uncertainties also tend to bereduced. The requirements of maximum parameter stability and minimum parameteruncertainty are balanced. The parameter estimates based on the set of weightssatisfying this criterion should be relatively stable, and have relatively small individualparameter uncertainties.1546.3.2.4 Procedure for weighting by the parameter space methods.The procedures for all three proposed parameter space methods are similar.The differences lie in what criterion is minimized during the selection of weights. Formaximum parameter stability, the CN is minimized. For minimum overall variance, thevolume of the parameter confidence region is minimized. The length of the longestaxis of the parameter confidence region is minimized as the third criterion.The basis for all of the parameter space methods is that each data setproduces an individual parameter confidence region. As the data sets are combinedusing different weights, the shape and orientation of the joint parameter confidenceellipsoid changes. For a particular set of weights, a confidence ellipsoid is obtainedthat minimizes the required criterion. The general procedure is:1. Minimize the objective function S(b) with w = w = 1.02. At minimum S(b), separate the approximate Hessian matrix into two parts,each part containing only head or only concentration information.H= wh Hh + w H (6.2)3. Scale the approximate Hessian matrices by the current parameter values.4. Scale the approximate Hessian matrices by __, where 52 is:s2 S(min) (6.3)n-pwhere S(min) is the minimum value of the combined objective function, n is the155total number of data, and p is the number of parameters being estimated. Thisscaling is done so that the inverse approximate Hessian matrix reflects thelinearized parameter confidence region. The s2 term also represents the overallvariance of the data.5. Varying the weights between 0 and 1, search for a set of weights thatminimizes the appropriate criterion. At each set of weights, decompose the totalHessian matrix into eigenvalues and eigenvectors that describe the parameterconfidence ellipsoid. The square roots of the inverses of the eigenvalues arethe lengths of the axes of the confidence region. These lengths are used tocalculate the criteria. The three criteria are:(1) MINCN: minimize the CN of confidence regionCN L(max) (6.4)L(min)(2) MINVOL: minimize the volume of confidence regionVOL = fl1L (6.5)(3) MINLEN: minimize the length of longest axisLEN = L(max) (6.5)6. Using the set of weights that minimized the appropriate criteria, estimate theparameters again. Repeat steps 2 through 6 until the set of weights stabilize.Less than three iterations are generally required for the weights to stabilize.156The proposed weighting methods produce unequal weighted residuals for eachdata set. This is not in accordance with standard regression techniques, whichdemand that the weighted residuals have equal variance. However, for joint data sets,the unequal variances of the weighted residuals may be necessary. Just as the dataand the model interact to define the parameter confidence region, the interaction ofthe data and model combine to determine the calculated variance of the dataresiduals. The variance of the data residuals for each data set are unique to that dataset, and this variance is somewhat model dependent. Due to this model dependence,unequal variances of the weighted residuals for each data type may produce the bestparameter estimates for a given model.Because the variances of the individual data sets are unequal, the assumptionsfor classical weighted regression are violated. The calculated parameter variancesmay not be equal to, and may underestimate, the true parameter variances. In thefollowing analysis, the parameter CV’S and uncertainties are calculated based onstandard weighted regression. The calculated CV’S are used to compare the differentweighting criteria. These calculated CV’s may not be equal to the true CV’S, and areused for comparison purposes only.6.3.3 Discussion of residual vs parameter space weighting criteriaWhat are the advantages of using the parameter space methods rather than157the analysis of data residuals method of selecting weights? The parameter spacemethods are designed to use more of the information available about the model beingcalibrated. The set of weights is designed based on the relative information each dataset carries about the model parameters. The set of weights reflects the informationeach data set has about all the parameters of the model.The parameter space methods incorporate information about the shape andorientation of the individual confidence regions in weighting the data sets. The dataresiduals approach determines the value of each data set (and therefore the weights)by the average residual variance for that data set. Two data sets may have similarresidual variances, yet their confidence regions may be quite different. If theconfidence region for head data is much larger than the confidence region forconcentration data, then the estimates using head data are probably further from thetrue values than the estimates using concentration data. The concentration data leadto better estimates, and should probably be weighted more heavily. In many cases,one set of data may provide better information on one of the parameters, while thesecond set of data provides better information on other parameters. The parameterspace methods allow the selection of weights based on the information each data sethas about all the model parameters.The parameter space methods can also adapt the weights to changes in themodel or parameter set being calibrated. If the model being calibrated is altered, or158the number of parameters being estimated is changed, the response surface for theparameter estimates will change. These different response surfaces result in differentconfidence regions, and parameter space methods can adapt the weights to theinformation contained in these confidence regions. In section 6.3.4.1, an example ispresented to show how the parameter space criteria adapt the weights to changes inthe model. In Chapter 7, a flow model is calibrated for two different parameter sets,illustrating how the parameter space criteria adapt the weights to different numbersof parameters being estimated.6.3.4 Comparison of weighting criteriaThe four weighting criteria can be compared using the synthetic flow andtransport system introduced earlier. Only two parameters are estimated; T1 and T2.The residual variances for each data set are approximately equal (h2 = 0.95; s2 =0.96). Gaussian random noise at a specified level of uncertainty was used to producethe data set, and the data covariance matrix is chosen to reflect these known datauncertainties. This choice results in residual variances for each data set which areapproximately equal to 1.0. Therefore, using the RESID criterion, the two data setswould be weighted equally. Figure 6.3a and 6.3b, introduced earlier, are the responsesurfaces for head and concentration data respectively.When the parameters are estimated with the joint data set, any of the four159-I 0 CD 0) F’)Q-13-4..-‘CD3DCDCD.c-‘S.CD(DCI)DCDCDD.- C) CD D CD a) -h 0 0 D-o 1:1) 3 CD -I. CD -S Co CD CI) Co -I’0 -SWeightingCriteria0) 0HeadDataTracerDataRESIDMINCNMINVOLMINLENWh1.000.001.001.000.651.00w0.001.001.000.211.000.40T1Estimate1.3031.3031.3031.3021.3031.302T2Estimate2.2842.2972.2972.2962.2972.297Std(T1)0.01020.01430.00790.00740.00820.0070Std(T2)0.04680.00890.00780.01210.00760.0097CV(T1)0.00790.01100.00630.00540.00630.0055CV(Ta)0.02040.00400.00340.00520.00320.0039AverageCV0.01420.00750.00490.00530.00480.0047ON2.943.011.871.242.071.39Volumeof2.OOe-044.30e-052.60e-052.95e-051 .99e-052.48e-05RegionLmax2.40e-021.14e026.90e-036.06e-036.44e-035.80e-03criteria for choosing weights may be used. Table 6.2 shows the weights, theparameter estimates and uncertainties for the four sets of weights. Figure 6.6a.,b.,c.and d. show the total response surfaces for the parameter estimates using the RESID,MINCN, MINVOL, and MINLEN criterion respectively.The MINCN criterion produces a parameter confidence region that is the mostcircular with the smallest CN, at the expense of larger parameter uncertainties (Figure6.6b). The parameter estimates based on the MINCN criterion have the largestvolume parameter confidence region, the longest length of the maximum axis, and thelargest average CV of any of the weighting schemes. The MINCN criterion weights thehead data set more heavily than the concentration data set in order to obtain the mostrounded confidence region. Because the head data set has larger parameteruncertainties the confidence region for the joint parameter estimates using the MINCNcriterion weights is large.The MINVOL criterion weights the concentration data set more than the headdata set (Figure 6.6c). This criterion results in a confidence region with the smallestvolume, but the parameter set has a larger condition number than any other criterion.The MINLEN criterion weights the head data set more than the concentration data set(Figure 6.6d), but not as severely as the MINCN criterion. The resulting parameterestimates have a relatively small CN and volume of the confidence region. Theaverage CV of the parameter estimates using the MINLEN criterion is the lowest. The161RESID criterion creates a parameter confidence region that is nearly as small as theMINVOL and MINLEN criterion. The parameter estimates for the MINLEN, MINVOLand RESID criterion are not significantly different.6.3.4.1 Adaptability to changes in the model -One way to illustrate the adaptability of the parameter space methods is tochange the model being calibrated. The same flow and transport system as above isused, and the same parameters (T1 and T2) are estimated. The model is changed sothat the inflow boundary is a specified flux boundary rather than a constant headboundary. The specified flux is calculated to produce the same head distributionthroughout the modelled area as the constant head boundary, and this flux value isspecified with no uncertainty. The data set used for parameter estimation is the sameas the original data set. This model is called the alternate model. The parameters T1and T2 are estimated, first using the head and concentration data sets individually. Theparameters are then estimated using the four weighting criteria to determine theweights for the joint data set.Figures 6.7a and b show the response surfaces for the individual data sets forthe parameter set T1-T2 using the alternate model. The parameter estimates and theconfidence regions are very different from the original model, even though the sameparameters are being estimated with the same data set. The data yield different162information about the model parameters because the model boundary conditions aredifferent. The head data alone produce much better parameter estimates than theconcentration data alone. In fact, the concentration data alone are virtually non-identifiable with respect to the two parameters being estimated. Because theconfidence regions are different, the set of weights required to minimize the parameterspace criteria are different. -Table 6.3 contains the parameter estimates, uncertainties, CV’S for theindividual data sets and the joint data sets using the four weighting criteria. Theresponse surfaces for the RESID, MINCN and MINVOL criteria are shown in Figure6.8a,b and c respectively. The RESID criterion yields the same set of weights (wh =1.0, w = 1.0) as the original model. However, the same set of weights does notproduce the same parameter estimates or confidence region for the joint data set. Thealternate model yields parameter estimates with smaller confidence regions for thejoint data set than the original model when using RESID or MINLEN weighting criteria.Using this alternate model, each of the parameter space criteria producedifferent sets of weights. The MINCN criterion (Figure 6.8b) weights the concentrationdata more than the head data, even though the concentration data alone yields a nonidentifiable parameter set. The shapes and orientations of the confidence regionsdictate that weighting the concentration data more than the head data results in a jointconfidence region that is as round as possible. The MINVOL criterion (Figure 6.8c)163p) CD 0) )(Do)D_DCD.CDCD-4. (Dp)Q-4-.CDo ciCDCDc 0 1 0 -4.ciCl) CD -4. C Cl) D Co -4’0WeightinqCriteria- 0)HeadDataTracerRESIDMINCNMINVOLMINLENDataWh1.000.001.000.051.000.98w0.001.001.001.000.101.00T1Estimate1.3031.2251.3011.2991.3031.301T2Estimate2.2992.2252.3002.3012.2992.300Std(T1)0.00867.07000.00590.00760.00590.0059Std(T2)0.00457.07000.00330.00650.00310.0033CV(T1)0.00665.90000.00450.00580.00450.0045CV(T2)0.00203.16000.00140.00280.00130.0014AverageCV0.00434.53000.00300.00430.00290.0030CN6.925884.582.216.324.58Volumeof6.60e-066.50e-014.65e-061.67e-053.66e-064.65e-06RegionLm6.70e-037.11e+004.60e-036.08e—034.87e-034.60e-03weights the head data more than the concentration data, resulting in a joint confidenceregion that has the minimum volume. However, the MINVOL criterion producesparameter estimates with the largest CN. For this alternate model, the MINLENcriterion results in the same set of weights as the RESID criterion. The MINLENcriterion and the RESID criterion result in parameter estimates that balance the needfor parameter stability and minimum volume of the joint parameter confidence region.This alternate model example simply demonstrates that even when estimatingthe same parameters with the same data set, the confidence regions can be verydifferent. Under the RESID weighting criterion, the weights are unchanged, but theconfidence regions are different. Under the parameter space criteria, the weightschange in response to the change in confidence regions. The weights changebecause the parameter space methods recognize the change in information availablefrom each data set.6.3.5 Discussion of results from tour weighting criteriaThe parameter space methods presented above offer an alternative to the dataresiduals weighting method for joint parameter estimation. The parameter spacemethods have the apparent advantage of incorporating information about parameterconfidence regions. From a conceptual standpoint, the MINLEN criterion would seemto result in the most sensible set of parameter estimates. It balances the need for165maximizing parameter stability and minimizing the total parameter uncertainty.In examples presented, the MINCN criterion often resulted in parameterestimates with the largest uncertainties. This was a direct result of the MINCN criterionminimizing the condition number at all costs. The MINVOL criterion did result inminimizing the volume of the joint parameter confidence region, but the parameterestimates often had the largest condition numbers. The MINLEN criterion generallyresulted in the parameter estimates with the lowest average coefficient of variation.The volume of the confidence region and the condition number of the parameter setwere both usually small.The RESID criterion also performed well, as measured by the averagecoefficient of variation, volume of joint confidence region, and condition number. TheRESID criterion seemed to do a good job in balancing the different criteria used toassess the parameter estimates. Based on the examples presented, as well as otherexamples not presented, the MINLEN and the RESID criterion seem to be the mostreasonable to use for weighting data sets in joint parameter estimation. Both criteriaresulted in parameter estimates which balance the need for maximizing parameterstability and minimizing parameter uncertainty.It may be appropriate to ask why the RESID criterion performs as well as theMINLEN criterion when it doesn’t include information from the parameter space. The166explanation may be that, for these lower dimensional parameter spaces, a wide rangeof weights produce a joint confidence region that changes in shape, but doesn’tchange the length of the longest axis significantly. As the weights change, the shapeof the joint confidence region does change, but this change is constrained by theshape and orientation of the individual confidence regions. Though the MINLENcriterion does minimize the length of the longest axis, this minimum length is not muchdifferent than the length for any reasonable set of weights. The RESID criterion resultsin a reasonable set of weights, and therefore does about as good a job as theMINLEN criterion. When more than two parameters are estimated, the differencesbetween the parameter estimates using the parameter space criteria and the RESIDcriteria should be greater. Examples of multi-parameter estimation are shown in thefollowing section and Chapter 7.6.4 Analysis of multi-parameter system using joint data setsA multi-parameter problem was examined using only hydraulic head data inChapter 5. The axes of the confidence region were analyzed to determine the mostefficient and responsible use of prior information to stabilize the parameter set. Ifconcentration data can be collected for this system, the additional data set may beable to significantly stabilize this parameter set. The contribution of a concentrationdata set is evaluated in this section.167The flow model and boundary conditions have been introduced in Chapter 5(Figure 6.14). To simulate a concentration data set, a constant concentrationconservative tracer is introduced along the western 2000 meters of the northernboundary of the recharge zones. For this simulation, all units are assigned an effectiveporosity of 0.2, and the specified transport parameter values are c=lOO m, T = 10m, and c = 1.0. This tracer is sampled at the 18 sampling locations at 5000 days.The standard error in the tracer concentration data is 5.0% of the concentration value.The concentration distribution within the flow system at 5000 days is shown in Figure6.9. The potential for this data set to reduce the non-identifiability of the 13 parameterset is to be evaluated.6.4.1 Analysis of parameter space based on concentration data setThe confidence region for the concentration data set needs to be defined. The13 flow parameters are estimated using only the concentration data. Using the scaledHessian matrix evaluated at the parameter estimates, the axes of the parameterconfidence region are given in Table 6.4. Table 6.4 shows that axis 13 of theconfidence ellipsoid is very long and nearly parallel to axis H. The orientation of thislongest axis of the confidence region indicates that H is an insensitive parameter. Themodel has no way of detecting the value or uncertainty of parameter H using theconcentration data. Table 6.5 lists the relative contribution to the uncertainty in eachparameter from each axis of the confidence region. Axis 13 contributes all the168uncertainty to parameter H, and contributes some of the uncertainty to the otherparameters. However, axis 12 also contributes significantly to the uncertainty of allparameters except H and R2.Table 6.6 contains the error ratios for parameters using the concentration dataset. The error ratios for the concentration data set are generally similar in magnitudeto the error ratios for the head data set, with the exception of errors in the specifiedhead boundary. Using concentration data only, errors in the specified head boundaryhave no influence on the estimates of the other parameters. However, errors in theother parameters may have a large influence on the estimate of the head boundaryvalue. The error ratios for the transmissivity parameters show no distinct pattern,unlike the error ratios using the head data only. The error ratios for the rechargeparameters, especially R1, are the largest, indicating that the estimates of the otherparameters are most sensitive to errors in recharge.In a multiparameter inverse problem, determining the difference in shape andorientation of the parameter confidence regions is more complicated than in the twoparameter problem. If only the longest axes were compared, axes 13 for both datasets would be used. However, for the concentration data, axis 13 contributes mainlyto the uncertainty in parameter H and gives little information on the uncertainty of theother parameters. Using head data, axis 13 gives information on the uncertainty of allother parameters except H. In this example, parameter H is weakly correlated with all1690>Dx0CDOu)DOCDDDO_4.-‘-‘ OCD -‘00-I-Q a- CD D 0 CD -S CD 0 D 0 -S 3 C 0 -, 0 3 CD CD -‘-o-S 0 0 CD 3 0 C,) CD 0-1 0 0 CD 0)Axis12e-021.7e-023.6e-024.7e-021.le-011.3e-011.6e-012.7e-014.5e-019.2e-011.2e÷001.le+012.2e+05LengthsParametersU1U2U3U4U5U6U7U8U9U10U11U12U13T1-0.300.21-0.020.72-0.080.49-0.110.26-0.050.10-0.070.100.00T20.010.070.160.01-0.68-0.44-0.140.26-0.040.340.160.310.00T30.06-0.07-0.18-0.170.540.07-0.280.230.180.420.260.480.00T4-0.400.26-0.230.180.29-0.60-0.36-0.08-0.25-0.060.00-0.200.00T50.16-0.040.43-0.21-0.010.28-0.630.05-0.48-0.080.08-0.160.00T60.18-0.31-0.260.040.09-0.090.270.69-0.41-0.010.02-0.260.00T7-0.510.20-0.39-0.59-0.260.32-0.020.13-0.040.07-0.08-0.090.00T8-0.210.360.63-0.150.27-0.090.250.360.220.17-0.09-0.240.00T90.14-0.21-0.140.10-0.090.02-0.26-0.090.330.57-0.18-0.600.00R1-0.59-0.750.260.020.01-0.05-0.02-0.020.07-0.050.010.060.00R2-0.100.020.080.020.080.080.39-0.41-0.550.560.190.000.00H0.000.000.000.000.000.000.000.000.000.000.000.001.00F0.07-0.050.03-0.060.08-0.08-0.030.00-0.200.15-0.900.320.000-1c5-CD CiicJ CDD (0..-I.DCD C:3_-..--O CDParametersU1U2U3U4U5U6U7U8U9U10U11U12U13T10.000.000.000.000.000.000.000.000.000.010.010.860.12T20.000.000.000.000.000.000.000.000.000.010.000.920.07T30.000.000.000.000.000.000.000.000.000.000.000.710.28T40.000.000.000.000.000.000.000.000.000.000.000.710.29T50.000.000.000.000.000.000.000.000.010.000.000.870.11T60.000.000.000.000.000.000.000.000.000.000.000.740.26T70.000.000.000.000.000.000.000.000.000.000.010.570.42T80.000.000.000.000.000.000.000.000.000.000.000.660.34T90.000.000.000.000.000.000.000.000.000.010.000.690.31A10.000.000.000.000.000.000.000.000.000.000.000.600.40R20.000.000.000.000.000.000.010.020.090.380.080.000.44H0.000.000.000.000.000.000.000.000.000.000.000.001.00F0.000.000.000.000.000.000.000.000.000.000.060.610.32-L-1 CD 0)-4.-5 0 3 0 -S >< -I. 0 3 C -4.-o 0 3 CD -4. CD -.5-o B CD 3 0• 0 Cl) CD 0 D C) 0 D C) CD D ‘-4. -5 0 ‘-4. 0ParameterswithPriorInformationEstimatedT1T2T3T4T5T6T7T8T9R1R2HFParametersT11.000.250.18-0.54-0.55-0.37-0.89-0.40-0.161.61-0.420.000.04T22.621.000.67-2.04-2.08-1.41-3.34-1.50-0.626.07-1.580.000.14T33.731.341.00-2.89-2.96-2.00-4.75-2.14-0.888.62-2.250.000.19T4-1.17-0.42-0.301.000.930.631.500.670.28-2.720.710.00-0.06T5-0.92-0.33-0.240.721.000.501.170.530.22-2.130.560.00-0.05T6-1.74-0.62-0.451.351.381.002.211.000.41-4.011.050.00-0.09T7-0.57-0.21-0.150.440.450.311.000.330.13-1.320.350.00-0.03T8-1.66-0.60-0.431.291.320.892.111.000.39-3.831.000.00-0.09T9-4.05-1.46-1.043.143.212.185.162.321.00-9.362.440.00-0.21R10.410.150.11-0.32-0.33-0.22-0.52-0.24-0.101.00-0.250.000.02R2-0.69-0.25-0.180.530.550.370.880.390.16-1.591.000.00-0.04H789.663853.913.1024.592455.828.933104411281.001241.F0.680.250.17-0.53-0.54-0.37-0.87-0.39-0.161.58-0.410.001.00rb,)other parameters for both data sets, but is well defined using head data alone andinsensitive using concentration data alone. The longest axes cannot be comparedblindly; a comparison of only the orientations of axes 13 for both data sets for thepurposes of determining if a reduction in parameter uncertainty will occur would be amistake in this case.Table 6.5 can be used to determine which axes of the confidence ellipsoidcontribute most to the uncertainty of the parameters. Axis 12 contributes the majorityof the uncertainty to all parameters except H and R2. Axis 13 contributes the majorityof the uncertainty to these two parameters. To determine if concentration data canincrease the stability of the inverse problem for the majority of the parameters, theangle between axis 13 (head data) and axes 12 and 13 for the concentration dataneed to be determined. Comparing axis 13(head) to axis 13(conc), 3 = 89.8°, andcomparing axis 1 3(head) to axis 1 2(conc), 3 = 82.3. The angle of interaction is quitelarge for both cases, indicating the concentration data may help stabilize and improvethe parameter estimates significantly.6.4.2 Joint parameter estimation using head and concentration data setsFor joint parameter estimation, weights must be selected for each data set.These weights can be selected on the basis of any of the four criteria given in section6.3. The weights for all criteria, along with the parameter uncertainties, CV’S, and CN’s173of the parameter sets, are given in Table 6.7. This table also gives the parameteruncertainties and CV’s for the head data set and the concentration data set alone.The average CV for all weighting criteria are given to allow easy comparison betweencriteria. The RESID criterion weights both data sets equally. The MINCN criterionweight the head data much less than the concentration data, while the MINLEN andMINVOL criteria weight the concentration data much less than the head data.The concentration data alone yields parameter uncertainties similar to thoseusing head data alone, except for parameter H, which has an extremely largeuncertainty. The large uncertainty of parameter H occurs because it is insensitive. Anycombination of data sets produce parameter estimates which are much better by anycriterion than the parameter estimates using individual data sets. The MINCN criterionresults in parameter estimates with a low condition number, but the CV of theparameter estimates is larger than any other criterion. This behavior is typical of theMINCN criterion, which maximizes parameter stability at the expense of increasing theuncertainty in the parameter estimates. The RESID criterion results in reasonablygood parameter estimates. However, both the MINVOL and the MINLEN criteriaproduce somewhat better parameter estimates, as measured by the averagecoefficient of variation. The MINLEN criterion results in the smallest average CV, yethas the largest CN of any of the criteria. It is difficult to single out one criterion as thebest, because the criterion which produce smaller average CV’S also produce thelarger CN’s.17403CD3-’(DC(DC)-‘CDo.TD30 -h 0 0 0 Cl) CD C C) D -h 0 C CD CD zr-I D CD C) CD -‘ -h 00 CD 0)-L-4 (51HeadDataOnlyConcentrationRESIDCriteriaMINCNCriteriaMINVOLCriteriaMINLENCriteriaDataWh1.00.01.00.021.01.0w00.01.01.01.00.190.10ParametersStdCVStdCVStdCVStdCVStdCVStdCVT15.762.597.152.490.110.050.450.220.090.040.080.04T25.662.7016.863.810.110.050.520.210.100.050.100.05T35.312.6512.312.890.230.120.700.300.180.090.170.09T45.532.422.411.220.110.050.210.100.080.040.070.03T55.463.242.561.480.050.030.120.070.040.030.040.02T65.284.393.732.720.040.030.090.080.030.030.030.03T75.523.074.672.290.060.040.180.100.040.020.040.02T85.334.415.793.740.020.020.080.060.020.020.020.02T95.257.155.997.560.010.020.040.060.010.020.010.02R13.8Oe-O12.707.Oe-050.2545.3e-060.01998.Oe-060.0251.00e-00.0251.00e-00.025R21.80e-O:12.631.6e-041.32.3e-050.1642.6e-050.2952.00e-00.1622.00e-00.165H1.750.021.8e÷061.8e+051.380.014.050.040.980.010.910.01F971.7713.73298.206.6216.400.2839.541.0011.540.1910.790.17Average5.5151.6e+050.0690.1950.0560.053CVConditionNumber1.le+041.8e+07134.74.144.148.VolumeofRegion1.32e-12.6e-0’2.17e-16.77e-15.43e-2(6.53e-20LengthofMaxAxis22.552.2e÷050.36951 .0020.30390.29926.4.3 Analysis of joint confidence regionUsing the scaled version of the approximate Hessian matrix containing bothhydraulic head and mass concentration data, evaluated at the parameter estimates,the axes of the joint confidence region can be calculated. Two aspects of this jointconfidence region relating to error ratios are evaluated; the influence of errors in thetransport parameters, and a comparison of the error ratios for the head and joint datasets.For this flow system, the concentration data reduce the uncertainty of theparameter estimates significantly. However, some assumptions were made whenconcentration data was simulated. Additional parameters were required to simulate theconcentration distribution. These parameter are the porosity of the internal zones, thesource strength, and the dispersivity of the medium. All of these parameters wereassigned a specified value without uncertainty. The effect of possible errors in theseparameter values on the parameter estimates needs to be determined. In particular,it would be nice to know if errors in these specified parameters will lead to large errorsin the estimated parameters. The linearized error ratio matrix for all parameters,including the transport parameters, must be calculated from the confidence regionusing the joint data set, even though these transport parameters only enter into thetracer data. If the linearized error ratios for any of the transport parameters are large,176then it may be better to estimate those parameters with large error ratios rather thanspecify them.The portion of the linearized error ratio matrix, calculated using equation (5.1),containing these transport parameters is included in Table6.8. The porosity of eachtransmissivity zone is defined as a separate parameter, (given as P1 through P9 inTable 6.8), but the longitudinal and transverse dispersivities are considered a singleparameter throughout the entire flow domain. It can be seen that errors in theconcentration source value contribute the most to errors in the estimated parameters.For the porosity parameters, the error ratios are generally small. However, the errorratios for lower transmissivity zones are larger than the error ratios for the highertransmissivity zones. In granular porous media, porosity is relatively easy to estimateindependently of the model, and large errors in estimates of porosity are uncommon.The dispersivity parameters have very small error ratios. Large errors in the specifiedvalue of dispersivity will only lead to small errors in the estimated values of the otherparameters. Estimates of dispersivities are scale dependent and usually very difficultto obtain independent of the model. The small error ratios for the dispersivityparameters means that poor estimates of the dispersivities will not introduce largeerrors in the estimates of the flow parameters. From this examination of the error ratiomatrix, the only transport parameter that may need to be estimated is the sourcestrength. We can cautiously state that the concentration data added significantly to theparameter stability and reduced parameter uncertainty without requiring accurate177-1 CD 9) cx m -‘ 0 0 3 -.‘ x 0 -S -S D-Cl)-o 0-D 3 CD CD -‘ Cl) C Cl) D CD 0 D -I 0 C’) CDCD Cl) 3 CD U) 0 - CD 0 0 a -I 0 D 0 -I -S 0 D Cl) 0 0 -S 0 3 CD CDTransportParameters-4EstimatedP1P2P3P4P5P6P7P8P9LTC1ParametersT10.00-0.02-0.010.010.020.06-0.010.060.090.000.00-3.8T20.00-0.01-0.010.000.010.040.000.040.060.000.00-4.0T30.00-0.02-0.020.010.020.07-0.010.070.100.000.00-3.5T40.00-0.01-0.010.000.010.040.000.040.060.000.00-5.8T50.00-0.02-0.010.010.020.06-0.010.060.090.000.00-6.8T60.00-0.02-0.020.010.020.08-0.010.080.120.000.00-4.7T70.00-0.02-0.010.010.020.06-0.010.060.090.000.00-6.9T80.00-0.02-0.020.010.020.08-0.010.080.120.000.00-11.1T90.00-0.04-0.030.010.040.14-0.020.130.210.000.00-58.2R1-0.01-0.07-0.050.020.070.23-0.030.220.36-0.01-0.01-18.0R2-0.01-0.07-0.060.020.070.24-0.00.20.3-0.01-0.0119.2H0.000.000.000.000.000.000.000.000.000.000.000.00F-0.01-0.09-0.070.030.090.30-0.030.290.47-0.01-0.01-18.9it is informative to compare the error ratio matrix for the parameter estimatesusing the joint data set to the error ratio matrix for the parameter estimates from thehead data set alone. The error ratio matrix for the joint data set is shown in Table 6.9.The error ratios are generally smaller for the joint data set than for the single data sets(Table 5.12 for head data; Table 6.6 for concentration data). The error ratios for theparameter F is the exception to this generalization. Parameter F has the mostuncertainty associated with its estimate during joint parameter estimation. The longestaxis of the joint parameter confidence region is nearly parallel to the axis of F inparameter space, resulting in large error ratios when estimating F. The parameterconfidence region for the joint data set is much smaller than the confidence region forthe single data set, and much better conditioned. These two factors allow the data setto compensate for errors in prior information. However, since this joint data set is well-conditioned for the 13 flow parameters, it is not imperative to use prior information toestimate the parameter. The two data sets together will do a good job in estimatingall the parameters. Prior information with uncertainty may be used to decrease thesize of the confidence region, if the uncertainty in the prior information is less than theuncertainty of the parameter from parameter estimation. Errors in prior information willnot have a significant effect on the final estimates, as these estimates are defined bythe well-conditioned response surface for the data.179H CD P CD ni -‘ -S 0 -S 0 3 a -S x 0 -S a CD 3 C ‘-I.-D-S 3 CD -4. CD -S-D -S 0 0 CD 3 C C’) D CD -4. CD 0 D a 0 (0 CD -4.ParameterswithPriorInformation-& 0EstimatedT1T2T3T4T5T6T7T8T9R1R2HFParametersT11.000.02-0.070.28-0.09-0.120.56-0.74-0.34-0.15-0.04-0.260.11T20.061.00-0.010.04-0.01-0.020.08-0.11-0.05-0.02-0.01-0.040.02T3-0.33-0.021.00-0.220.070.10-0.450.600.270.120.030.21-0.09T40.240.01-0.041.00-0.05-0.070.33-0.44-0.20-0.09-0.02-0.150.06T5-0.050.000.01-0.031.000.01-0.070.090.040.020.000.03-0.01T6-0.060.000.01-0.040.011.00-0.080.110.050.020.010.04-0.02T70.290.01-0.050.19-0.06-0.081.00-0.52-0.24-0.10-0.03-0.180.07T8-0.14-0.010.02-0.090.030.04-0.181.000.110.050.010.09-0.03T9-0.050.000.01-0.030.010.01-0.070.091.000.020.000.03-0.01R1-0.040.000.01-0.020.010.01-0.050.060.031.000.000.02-0.01R2-0.63-0.030.10-0.420.140.18-0.851.140.510.221.000.40-0.16H-0.030.000.00-0.020.010.01-0.040.050.020.010.001.00-0.01F3.860.19-0.612.58-0.83-1.105.21-6.97-3.15-1.36-0.34-2.451.006.5 SummaryThree aspects of joint parameter estimation have been discussed in thischapter. First, response surfaces and confidence regions were used to show howmultiple data sets reduced parameter uncertainty. Second, the use of confidenceregions was extended to predict the value of future data in reducing the uncertaintyof parameter estimates, before the data are collected. Third, parameter spaceapproaches were introduced for selecting the weights for the individual data sets injoint parameter estimation. Four weighting criteria were compared using both simpletwo parameter problems and a multi-parameter problem.Each data set produces a unique response surface and confidence region. Ifthese surfaces are oriented differently for the two data sets, the parameter uncertaintyare significantly reduced when the second data set is included. If these surfaces aresimilar for two data sets, the second data set will not reduce the parameter uncertaintysignificantly. In simple two parameter examples, the response surfaces werevisualized. For multiparameter examples, the axes of the confidence ellipsoid need tobe analyzed to determine the difference in orientation of the confidence regions for thetwo data sets.Future data can be simulated using the model and the parameter estimatesfrom the original data set. The confidence region for the future data set can be181estimated using these simulated data values. The confidence regions forthe initial andfuture data sets can be compared to determine the potential reduction in parameteruncertainty for the future data set. Though the actual future data will not produce thesame parameter estimates as the simulated future data, the actual future data willproduce similar values for parameter uncertainty.For weighting data sets in joint parameter estimation, the analysis of dataresiduals method was compared to three methods based on parameter spaceanalysis. The parameter space based methods have the advantage of incorporatinginformation about parameter confidence regions. Three criteria based on parameterspace analysis have been proposed (minimum condition number, minimum volume ofconfidence region, and minimum longest axis of confidence region). Each criterionresults in the “best” parameter estimates in some sense. The parameter space criteriacan adapt the weights to changes in the model. The MINCN criterion often resultedin parameter estimates with the largest uncertainties, but lowest condition number.The MINVOL criterion resulted in parameter estimates with the smallest sized jointconfidence region, but the estimated parameter set often had larger condition numberscompared to other weighting criteria. In the examples presented, the MINLEN criterionand the RESID criterion seemed to be the most reasonable to use for weighting datasets in joint parameter estimation. They both balanced the need for maximizingparameter stability and minimizing the total parameter uncertainty.182A multi-parameter flow system was used to demonstrate the reduction inparameter uncertainty during joint parameter estimation. The parameter set was non-identifiable using only head data. Tracer concentration data significantly reducedparameter uncertainty and ill-conditioning. Additional transport parameters wereneeded to simulate the tracer concentrations, but errors in the specified values ofthese parameters were found to lead to very small errors in the estimates of the flowparameters, based on linearized error ratio analysis. The concentration data addedsignificantly to the parameter stability and reduced parameter uncertainty withoutrequiring accurate estimates of the transport parameters. With respect to weightingthe data sets, it was difficult determine which criterion performed the best. TheMINLEN criterion resulted in the smallest average coefficient of variation, but it alsoresulted in the largest condition number. The RESID criterion, MINVOL criterion andMINLEN criterion all resulted in reasonably good parameter estimates.18311 Co C -.‘ CD‘<Q)0,0CD -‘ 0 0 D 0 CD -‘ CD 0 :3 0 C,,ODc4-. 0 :3 0 Cli-o cx) C,) CD CD 0 (I) C)) CD C,) (ii 0 0Sourcearea—I6000m1.151.091.03R 0.970.910.850.85 0.91 0.97 1.03 1.09 1.15(a)1.151.091.03R 0.970.910.850.85 0.91 0.97T1(b)1.03 1.09 1.151.151 .091.03R0.970.910.850.85 0.91 0.97 1.03 1.09 1.15(C)Figure 6.2 Response surfaces for parameter set A: (a) head data; (b) concentrationdata; (c) joint data set1851 .051 .031.01T2 0.990.970.951.051 .031.01T2og0.970.951.051 .031.010.990.970.950.95 0.97 0.99 1.01T1(C)1.03 1.05Figure 6.3 Response surfaces for parameter set B: (a) head data; (b) concentrationdata; (c) joint data set0.95 0.97 0.99 1.01 1.03 1.05 0.95 0.97 0.99 1.01 1.03 1.05T1(a) (b)I I I I I I I1861.05 1.051.03 — 1.03 —1.01 1.01 —T20.99 T2 o.0.97 0.97 —0.95 I I I II c.95 I I I I I I I0.95 0.97 0.99 1.01 1.03 1.05 0.95 0.97 0.99 1.01 1.03 1.05T1(a) (b)1.05 1.051.03 1.03 -1.01-1.01 -::: :0.95 I I I I I I 0.95 I I I I I0.95 0.97 0.99 1.01 1.03 1.05 0.95 0.97 0.99 1.01 1.03 1.05T1 T1(C) (d)Figure 6.4 Response surfaces for T1-T2 parameter set: (a) synthetic futureconcentration data; (b) head and synthetic future concentration data; (c)actual concentration data; (d) head and actual concentration data187(a)(b)Figure 6.5 Two confidence ellipses with the same volume and orientation, butdifferent stabilities and parameter uncertainties1881 .05 1 .051.03 -0.950.95 0.97 0.99 1.01 1.03 1.05T1(a)1.031.01T2•g0.970.950.95 0.97 0.99 1.01(b)I I I I I I I I1.03 1.050.970.950.95 0.97 0.99 1.01 1.03 1.05Figure 6.6 Response surfaces for T1-T2 parameter set using joint data sets: (a)RESID weights; (b) MINCN weights; (c) MINVOL weights; (d) MINLEN1.01 —0.97 —I I I I I I I IT11.03 1.051.051.031.010.991.051.03 -1.010.990.970.950.95 0.97 0.99 1.01(C)I I I I I I(d)weights 1891 .051.031.01T20.990.970.950.95 0.97 0.99 1.01T1(a)1.03 1.051.051.031.01T2 0990.970.950.95 0.97 0.99 1.01 1.03 1.05(b)Figure 6.7 Response surfaces for alternate model: (a) head data; (b) concentrationdataT11901.05 1 .050.97—0.95 —0.95 0.97 0.99 1.01 1.03 1.050.970.950.95 0.97 0.991.051.031.010.990.970.950.95 0.97 0.99 1.01 1.03 1.05(C)Figure 6.8 Response surfaces for alternate model using joint data sets: (a) RESIDweights; (b) MINCN weights; (C) MINVOL weights1.03 —T20991.031.010.99(a)I I I I I I I I1.01 1.03 1.05(b)ENI I I I I191600050004000300020001 0000 6000Figure 6.9 Concentration distribution for multi-parameter flow system at 5000 days0 1000 2000 3000 4000 5000192CHAPTER 7. JOINT PARAMETER ESTIMATION FOR THE SAN JUAN BASINUSING PARAMETER SPACE METHODSIn this chapter, the development and calibration of a groundwater flow modelfor the San Juan Basin, New Mexico, is discussed. Hydraulic head data, 14C data, andprior information on parameter values are used for parameter estimation in the SanJuan Basin. A cross-sectional model of the hydrogeological units in the west-centralportion of the basin is constructed. The purpose of constructing the model is todemonstrate the utility of the methods developed in this thesis for calibrating ahydrogeological model in an efficient and responsible manner using joint data sets.The geology of the basin has been investigated by Baltz [1967] and Fasset and Hinds[1971]. The hydrogeology of the basin has been summarized by Stone et al. [1983].7.1 Study areaThe San Juan basin is a large structural basin (77000 km2) in northwesternNew Mexico and southwestern Colorado (Figure 7.1). The basin lies in a relatively flatupland region south and east of the Rocky Mountains. Altitudes range from 2400 min the northern part of the basin to about 1550 m where the San Juan river exits thebasin. The central part of the basin is a dissected plateau, the surface of which slopesgently to the west. Figure 7.1 is a plan view of the San Juan Basin, showing thesurface outcrop of the major Cretaceous and Tertiary units. The line A-A’ in the193southwestern portion of the basin is the line of section along which the model isconstructed. The model extends from Chaco Canyon in the south to the San Juanriver in the north, along the center of R1OW.The climate throughout-the basin is arid to semi-arid. The rainfall in the basinis about 10 cm a year at the lower altitudes of the modelled cross section. The majorstream in the basin is the San Juan River, which flows westward in the New Mexicoportion of the basin. Many intermittent streams also drain the basin. The Chaco Riveris the largest intermittent stream, draining the southern and western parts of the basinand flowing into the San Juan on the western side of the basin.7.2 GeologyThe San Juan basin is nearly circular in plan view. It is an asymmetric structuraldepression which contains sedimentary rocks ranging in age from Cambrian toQuaternary. This study will concentrate on groundwater flow in the upper Cretaceousand Tertiary units of the basin. Figure 7.2 is a schematic cross section of the upperCretaceous and Tertiary units along the modelled section.The upper Cretaceous rocks of the basin are over 1830 m thick and consist ofintertonguing marine and non-marine units that represent three basin-widetransgressive-regressive cycles of deposition. These cycles laid down the Menefee194Formation, the Cliff House Sandstone, the Lewis Shale, with the final regression of thesea resulting in the deposition of the Pictured Cliffs Sandstone. The Pictured Cliffs isoverlain by two Cretaceous units of continental origin, the Fruitland Formation and theKirtland Shale. The Tertiary rocks in the study area consist of the Ojo AlamoSandstone and the Nacimiento Formation. The Nacimiento contains interbeddedmudstone with sandstone lenses. Quaternary deposits consist of Pleistocene andHolocene terrace gravels and alluvium.7.3 HydrogeologyThree types of hydrogeological information are available to construct ahydrogeologic model of the basin; hydraulic head data, 14C data, and direct informationon parameter values.The head data were obtained from Stone et a!. [1983]. Phillips et a!. [1989]kriged 24 head data in the Ojo Alamo and 26 head data in the Nacimiento to calculatethe potentiometric surface for each aquifer. The head distribution indicates rechargealong the outcrops and discharge to the San Juan River. The head distribution alsoshows that the hydraulic head contours are generally perpendicular to the modelledcross section.For this study, all head data available in R9W, Ri OW and Ri 1 W from Chaco195Canyon in the south to the San Juan river in the north were analyzed for use in thecross section model. The sample locations are plotted on Figure 7.3, along with theinformation about which formation each well sampled. Most of the wells werecompleted in the Ojo Alamo and Nacimiento formations, but several were completedand sampled in lower formations. The screened intervals for the wells range from 5mto over 1 OOm, with an average of approximately 20m. Table 7.1 contains the well data -used in this study, listing a well identifier, well location, the elevation of the midpointof the screened interval, and the hydraulic head. The well identifier consists of anumber and letters. The letters indicate which formation the well was reported to havesampled. The letter identifiers are:(q) Quaternary Alluvium(n) Nacimiento(0) Ojo Alamo(k) Kirtland/Fruitland(pc) Pictured Cliffs(I) Lewis(ch) Cliff HouseThe hydraulic head data set consists of all wells in R9W, R1OW, and R11W, up to 15km from the line of section. Figure 7.4 shows the modelled cross section, along withthe midpoint of the screened interval for all wells in the head data set.The assigned uncertainty in the head data is determined in part by the distance196Table 7.1 Hydraulic headmodeldata used for calibrating the San Juan Basin35n39n48n51n52n27443.0 1976.0 1986.0 30.0Well Distance from Midpoint of Hydrauhc StandardIdentifier Chaco Canyon Screened Head Error(m) Interval (m) (m)(m)29n 30323.0 1993.0 2018.0 35.042181.0 1833.0 1953.0 30.040n 42011.0 1860.0 1927.0 25.046n 54547.0 1857.0 1855.0 45.047771.0 1833.0 1908.0 55.052683.0 1822.0 1838.0 55.047093.0 1775.0 1855.0 50.056n 62847.0 1840.0 1840.0 25.059n 67252.0 1738.0 1747.0 40.062n 69115.0 1654.0 1669.0 60.072n 75000.0 1658.0 1660.0 30.073n 74875.0 1652.0 1678.0 30.0190 25749.0 1906.0 2026.0 60.0200 25241.0 1915.0 2009.0 60.0210 25749.0 1909.0 1920.0 60.0300 28629.0 1928.0 1932.0 45.0310 33541.0 1960.0 1975.0 20.0320 33711.0 1857.0 1912.0 30.01856.0 1893.0360 33711.0 1920.0 1920.01845.0 1930.0380 40148.0 1810.0 1887.0410 40487.033034037°4n32017.027443.040656.07ngq n1924.01870.011Afl1930.01935.01R9fl20.020.055.035.035.020.020 0197Table 7.1 (continued)Well Distance from Midpoint of Hydraulic StandardIdentifier Chaco Canyon Screened Head Error(m) Interval (m) (m)(m)440 39809.0 1846.0 1876.0 55.0450 45230.0 1812.0 1850.0 60.0520 47093.0 1773.0 1855.0 50.0610 69285.0 1625.0 1658.0 55.0740 74028.0 1526.0 1654.0 60.053g 46416.0 1880.0 1880.0 45.068g 74738.0 1663.0 1663.0 30.069g 74399.0 1679.0 1683.0 30.070g 74044.0 1710.0 1719.0 20.04g 4743.0 1949.0 1946.0 30.016g 8131.0 1874.0 1874.0 30.017g 7623.0 1871.0 1871.0 30.025g 17448.0 1957.0 1956.0 30.026g 17787.0 1957.0 1956.0 30.08k 11180.0 1862.0 1862.0 50.010k 15754.0 1839.0 1893.0 20.012k 14568.0 1902.0 1915.0 20.013k 12197.0 - 1876.0 1894.0 25.014k 12366.0 1872.0 1891.0 30.015k 11858.0 1833.0 1880.0 30.020k 25241.0 1825.0 1935.0 60.027k 25579.0 1320.0 1787.0 60.058k 57765.0 1576.0 1737.0 40.0ich 4913.0 1826.0 1898.0 35.02ch 4913.0 1664.0 1852.0 50.06pc 14399.0 1815.0 1938.0 60.09ch 11350.0 1619.0 1731.0 50.0llpc 14399.0 1798.0 1882.0 25.0lRch 11180.0 18210 1849.0 45.0198of the well from the line of cross-section. This uncertainty has several sources, all ofwhich contribute to the modelled data riot matching the measured data exactly. Someof these sources are: (1) an uncertainty about which vertical point in the aquifer issampled due to the length of the screened interval of the well; (2) an uncertainty dueto the fact that-the wells do not lie on the line of section; and (3) the uncertainty dueto the natural variation in hydraulic conductivity and other parameters within each unit,which is present in the field but not included in the model. For wells within 0.8 km (onesection) of the line of the cross-section, the standard error of measurement is 20meters. For wells further from the line of the cross-section, the standard error ofmeasurement is 20 meters, plus 5 meters for each section (approximately 1.6 km)away from the cross section. This method of assigning uncertainty based on thedistance from the modelled section is rather arbitrary, but qualitatively seems to be areasonable method of assigning uncertainty in a cross sectional model. The baselineuncertainty of 20 m may seem to be large for measurements of hydraulic head, butit is consistent with the large scale of the model. Table 7.1 also contains the assignedstandard error of measurement for each well.The 14C data were obtained from Phillips et al. [1989] and Shute [personalcommunication, 1994]. These data are available only in the Ojo Alamo and Nacimientoformations. The 14C data available in R9W, R1OW, and Ri 1W were used for thisstudy. Table 7.2 contains the ‘4C data, reported in percent modern carbon, along withthe well locations and midpoint of the sampling interval. Some of the wells sampled199Table 7.2 The 14C data used to calibrate the San Juan Basin modelWell Shute Phillip Distance Midpoint C-i 4 c-i 4 [Phillip] StandardIdentifier Well Well from of [Shute] (pmc) ErrorNumber number Chaco Screened (pmc)Canyon Interval(m) (m)28n SJB-24-O 35405 1934 51.2 ± 0.6 9.139n NMO5 SJB-22-N 42181 1833 12.6 12.3 ± 0.5 2.340n NM22 42011 1860 19.9 2.548n NMO7 SJB-35-N 47771 1833 8.39 9.2 ± 0.4 4.481n NMO1 SJB-21-N 67500 1553 1.4 <1.73 0.582n NMO8 57700 1681 4.49 3.9200 NM1O 25241 1915 37.9 7.822o SJB-13-O 25750 1909 iO.16±.56 5.0230 NM17 SJB-04-O 17800 1934 50.5 51 .22± 1.14 7.5320 NMO9 SJB-02-O 33711 1857 27.2 28.66 ± .85 3.7330 SJB-17-O 32017 1856 40.3 ± 0.9 4.0340 SJB-05-O 27443 1924 5.52 ± .85 0.5380 NM11 SJB-08-O 40150 1810 8.33 8.99 ± .9 2.3520 NM13 47093 1773 8.14 3.8570 NMO3 62800 1626 2.79 3.7580 NMO2 SJB-15-O 57765 1624 3.96 10.25 ± .75 2.4830 NM23 SJB-12-O 51500 1700 8.85 4.40 ± 1.5 1.8840 SJB-11-O 47000 1756 4.70 ±0.6 1.4by Phillips et al. were re-sampled by Shute, and the ‘4C concentrations are inagreement for most of the wells. Two wells sampled in both studies have different ‘4Cconcentrations: wells 580 and 830. In this study, the more recent sampleconcentrations are used for model calibration.200Figure 7.5 is a location map of the wells from which the 14C data was collected.Eight wells lie within R1OW, and the remaining 10 wells are further from the line ofsection. Figure 7.6 shows the modelled cross section, along with the locations of themidpoint of the sampling interval for the wells. Some of the same wells sampled forhead data are also sampled for ‘4C data. For most of the wells, the midpoint of thescreened interval for hydraulic head data is identical to the midpoint of the samplinginterval for the concentration data. Two wells are different. Well 340 (SJB-05-O) isreportedly screened in the Ojo Alamo, but Phillips eta!. [1989] stated that the waterchemistry indicated that it came from the Kirtland Shale. For the model calibration, thesampling interval is located in the Kirtland Shale. Well 58 (SJB-15-O) was screenedas deep as the Pictured Cliffs, but was probably sampled in the Ojo Alamo, so thesampling interval is placed in the Ojo Alamo.The assigned uncertainty in the 14C data is determined using two factors: (1)the magnitude of the 14C concentration and (2) the distance of the well from the lineof cross-section. The sources of uncertainty in the 14C data are the same as those forthe hydraulic head data, as well as the measurement uncertainty reported by Phillipset al. [1989]. It is important to link the uncertainty in the 14C concentrations to themagnitude of the concentrations in order to capture the pattern of concentrationdistribution, If the uncertainty were not linked to the magnitude of the concentrations,the match between the simulated concentrations and the observed concentrationswould be dominated by the high concentration samples. The lower concentration201samples need to be matched more closely than the higher concentration samples inorder to capture the overall pattern of the concentration distribution. The assigneduncertainty for each sample is calculated as 10% of the ‘4C concentration, plus 0.5pmc for each section (approximately 1.6 km) away from the cross section. Table 7.2also contains the assigned standard error of measurement for each well.The direct information on hydraulic conductivities is from Stone eta!. [1983].This direct information can be used as either prior information on parameter valuesor initial parameter estimates for model calibration purposes. However, the priorinformation may not be representative of the model parameter values for severalreasons. First, the prior information is obtained from drill stem tests, slug tests, andpump tests, which sample the subsurface at a much smaller scale than the modelparameter zones. The prior information may underestimate the model parametervalues due to differences in scale. Second, the prior information in low permeabilityunits is often obtained in the more permeable sandstone lenses within a lowerpermeability mudstone or shale unit. The sandstone lenses are often notinterconnected, and the prior information obtained from the lenses may overestimatethe overall permeability for the units. Third, the prior information is often from locationsfar from the modelled area. This prior information is often taken near the outcrop, andmay not be representative of parameter values at depth. Table 7.3 contains the priorinformation available on the hydraulic conductivity values for the units in the model,along with the source of the prior information and the estimation method. Several prior202values are available for some of the units.A second source of prior information is the model of Phillips and Tansey[1 984].They constructed a quasi three-dimensional groundwater flow model of the Ojo Alamoand Nacimiento formations, using 14C ages to obtain the hydraulic conductivitydistribution. They divided the Nacimiento into three units; an aquifer and two aquitards.The thickness of the aquifer was obtained by adding the thicknesses of the sandstonelenses throughout the Nacimiento together, and considering this total thickness anaquifer located in the center of the Nacimiento. An aquitard was located above andbelow this aquifer. The model was calibrated by adjusting the permeability of theaquitards. Table 7.3 also contains the prior information available from the Phillips andTansey [1984].7.4 Model constructionA cross-sectional model of the hydrogeological units in the west-central portionof the basin, representing the southern limb of the basin, is constructed. The purposeof constructing the model is to demonstrate the utility of the methods developed in thisthesis for calibrating a hydrogeologic model in a responsible manner. Hydraulic headcontours for the basin show that the flow direction is generally N-S, so the model isconstructed along a north-south cross section through the middle of R1OW. Thesouthern boundary of the model coincides with the bottom of Chaco Canyon, and the203Table 7.3 Prior information on the parameters in the San Juan Basin modelFormation Prior Hydraulic Source Test MethodConductivityNacirniento 0.63 rn/yr Phillips and thermal gradientsaquitard Tansey [1984]Nacimiento 0.015 rn/yr Phillips and model calibrationaquitard Tansey [1984]Nacimiento 6-15 rn/yr Phillips and calculated from flowaquifer Tansey [1984] path and 14C dataNacirniento 11 rn/yr Stone et al., estimatedaquifer [1983]Ojo Alamo 15-50 rn/yr Phillips and calculated frorn flowTansey [1984] path and 14C dataOjo Alamo 55-177 rn/yr Brimhall [1973] aquifer pump testsFruitland/Kirtland 0.001-0.01 rn/yr Stone et al. [1983] estimated based onvariety of testsFruitland/Kirtland 0.015 rn/yr Phillips and thermal gradientsTansey_[1984]Pictured Cliffs .001-3.3 rn/yr Stone et al. [1983] aquifer testsPictured Cliffs 0.7 rn/yr Stone et al. [1983] permeability data fromReneau and Harris[1953]Lewis Shale - - no dataCliff House 1.11 rn/yr Stone et al. [1983] recovery testCliff House 0.11 rn/yr Stone et al. [1983] permeability data frornReneau and Harris[1953]Menefee 1 .Oe-6 to 1 .1 Stone et al. [1983] aquifer testsm/ynorthern boundary is located at the San Juan River. The distance from Chaco Canyonto the San Juan River is approximately 75 km along R1OW.2047.4.1 Model boundary conditionsThe locations of the model boundaries were chosen so that simple boundaryconditions could be applied. Figure 7.7 is a conceptualized cross section showing thelayering and boundary conditions. The flow system is assumed to be at steady state.1. Upper BoundaryThe upper boundary of the model is chosen to coincide with the elevation ofthe water table in the uppermost unit. From south to north along the cross section, thesurface elevation starts at 1820 meters at the base of Chaco Canyon, and increasesgradually to the outcrop of the Ojo Alamo, at approximately 2000 m. A plateaucontinues for about 15 km north of the outcrop of the Ojo Alamo, and then surfaceelevations decrease to the San Juan River, at approximately 1690 meters.The upper boundary is represented as a specified head boundary. To imposethis boundary, four nodes are located along the top of the flow system, one at ChacoCanyon, one at the outcrop of the Ojo Alamo (20 km), one at 35 km, and the fourthat the San Juan River (75 km). The value of the hydraulic head at each of thesenodes is estimated, and the hydraulic head along the rest of the boundary isinterpolated linearly between these nodes. During model calibration, the elevations ofthe top boundary change, and the FEM grid can adjust to reflect these changes. Theinitial estimate of the head along the top boundary is the surface elevation.2052. Lower BoundaryThe lower boundary is represented by a no-flow boundary within the MenefeeFormation. The Menefee Formation is the unit below the Cliff House, and is arelatively thick sequence of marine shales and siltstone lenses. Though morepermeable formations do exist below the Menefee Formation, the flow across theMenefee represents a small percentae of the total flow through the system. Thehydraulic conductivity of the Menefee is assigned a value of 1 .Oe-6 m/yr, and is notestimated. Several factors were taken into account in order to locate the domainboundary below the Cliff House. First, no hydraulic head or 14C data was available inunits below the Cliff House near the cross section, so the parameters for lower unitswould have been extremely difficult to estimate. Second, numerical experimentsshowed that including the units below the Cliff House in the model resulted in verylittle change in the hydraulic head or 14C concentration distribution in the upper layers.Including the lower units was not necessary in order to simulate flow and transport inthe units where data exists.3. Southern Boundary.The southern boundary of the model is chosen to coincide with the bottom ofChaco Canyon. This canyon is a suriace water and shallow groundwater divide. It maynot be a groundwater divide for the deeper regional groundwater system. However,the units intersecting this boundary are low permeability, so the no-flow boundaryspecification is probably reasonable.2064. Northern BoundaryThe northern boundary is specified to be a no-flow boundary. The San JuanRiver marks the lowest elevation (in a N-S cross section) of the basin, and also marksthe axis of the basin. The San Juan River should be both a shallow groundwater anda regional groundwater divide, so a no-flow boundary is appropriate.7.4.2 Hydrogeological units in the modelSix major hydrogeological units are included in the model. The vertical locationand dip of these units were obtained by plotting all data from the well logs, and fittingthe units to the well log completion intervals. The distribution and thickness of thehydrogeological units in the model is shown in Figure 7.7. The Cliff House Sandstone,the lowest stratigraphic unit estimated using the model, crops out in Chaco Canyonat an elevation of approximately 1820 m. Based on well data, it dips to the north atan angle of about 1.50. It is represented as a 50m thick unit. All other units dip to thenorth at an angle of about 10. The Lewis Shale crops out just north of Chaco Canyon.It is relatively thin near the outcrop and increases in thickness northward up to 525m.The Pictured Cliffs Sandstone crops out about 8km north of Chaco Canyon, andmaintains a constant thickness of 25m throughout the model domain. TheKirtland/Fruitland Shale is represented as a unit 200 meters thick. The Ojo AlamoSandstone crops out about 15 km north of Chaco Canyon, and dips north at 10. It isrepresented as a 50 m thick layer. The Nacimiento Formation is the highest207stratigraphic unit in the cross section. It is present only north of the outcrop of the OjoAlamo, and ranges in thickness from 0 to 150 meters. For this model, the Nacimientois considered a single hydrogeologic unit. Though the Nacimiento containsinterbedded mudstone and sandstone lenses, it is not known whether the sandstonelenses are interconnected. Defining the Nacimiento as a single hydrogeologic unit ismore reasonable than breaking it into several units, since the location and extent ofmultiple units would be difficult to determine.7.4.3 Concentration boundary conditions and parametersIn order to simulate the 14C data, the boundary conditions and parameters thatare specific to the tracer concentration field are required. The ‘4C is introduced alongthe top boundary of the model, where it enters the groundwater flow system at thewater table. This boundary is a specified concentration boundary. The requiredparameters are the half life of the 14C isotope, the initial concentration at the inflowboundary, the horizontal and transverse dispersivities, and the effective porosity foreach unit.The determination of the initial concentration of 14C at the water table is notstraightforward. Phillips et al. [1989] used six different correction models to estimatethe initial concentration of 14C samples in the San Juan basin, and obtained valuesranging from 31 to 150 prnc. Phillips eta!. [1989] determined that the correction model208of Vogel [1967] compared favorably with the other correction models they examined.The correction model of Vogel assigns an initial activity of 85 pmc to every sample.For the model in this thesis, an initial concentration of 85 pmc is used at the watertable. This boundary condition was applied along the water table in the region wheredownward flow from the water table to the aquifer was present. In the dischargeportion of the upper boundary, a third-type boundary was applied.The horizontal and transverse dispersivities are not straightforward todetermine. Phillips et al. [1989] used well logs and stochastic transport theory toestimate horizontal dispersivities of approximately 300 meters. However, using a onedimensional transport model, they showed that the concentration distribution was notvery sensitive to the horizontal dispersivity for a decaying solute. Comparing horizontaldispersivities ranging from 0 to 1500 m, the maximum difference in 14C ages wasfound to be 11%, within the measured uncertainty in the 14C concentrations. For themodel in this thesis, a horizontal dispersivity of 300 meters is used. The transversedispersivity is estimated at 10 meters. The effective porosity is assigned a value of0.2. After the parameters are estimated using the joint data set, the consequences ofincorrectly specifying the parameters controlling the concentration distribution areexamined.The half life of 14C is approximately 5740 years [Fontes, 1980]. Since the flowsystem is at steady state and is presumed to have been relatively unchanged for a209long period of time relative to the half life of 14C, the 14C concentrations are also atsteady state. The time needed to reach a steady state distribution depends on thedistance from the source. For a wide range of aquifer properties, it was determinedthat concentrations were relatively unchanged within the model domain after 50,000years. The maximum simulation time is specified to be 50,000 years.The ‘4C distribution is calculated both using standard Crank-Nicolson timestepping and the Arnoldi algorithm. Using the Crank-Nicolson time stepping procedure,relatively accurate steady state distributions can be calculated using time steps of 500years, so only 100 time steps are needed to calculate the 14C distribution. Using theArnoldi algorithm, the amount of computer time required to accurately calculate thesteady state 14C distributions is dependent on the contrast in hydraulic conductivitybetween the model layers. In general, at least 10 Arnoldi vectors are needed formoderate contrasts, and up to 15 Arnoldi vectors are needed when the contrast inhydraulic conductivity between layers is greater than 3 orders of magnitude. Inaddition, many more iterations are needed to solve for each Arnoldi vector when thecontrast between the layers is large. Because the hydraulic conductivity of each layeris variable during the parameter estimation procedure, and contrasts of greater thansix orders of magnitude are common, the Arnoldi algorithm performs poorly duringparameter estimation for this model. The Crank-Nicolson procedure is more robust,and hence more reliable during parameter estimation. It is often faster as well.2107.4.4 Model gridThe model is based on a finite element technique using triangular elements.The mesh is aligned with the regional slope of the hydrogeological units. Thisalignment allows easy representation of the units, and allows the elements to bealigned with the streamlines in the more permeable units. The upper boundary of themesh is adjustable so that it can coincide with the elevation of the water table in theupper units. For the flow model, the elements are 250 meters long and between 5 and15 meters high.7.5 Model calibration strategyThe flow model, as described above, contains 10 parameters that need to beestimated. There are 4 head boundary parameters; the value of hydraulic head ateach of the 4 nodes that define the upper surface of the cross section. There are also6 hydraulic conductivity parameters, one for each layer in the model. The anisotropyfor each model hydraulic conductivity zone also needs to be defined, but thisparameter is not estimated initially. The anisotropy is assigned a value of 100 for eachmodel layer, and subsequently the effect of errors in the specified value of anisotropyon the parameter estimates are determined.The model is calibrated in two stages. In the first stage, all 10 model211parameters are estimated. The hydraulic head and 14C data are used both individuallyand jointly to estimate the model parameters. The parameter space for the joint dataset is examined to determine how to use prior information to stabilize the parameterset most efficiently and responsibly. It is shown that the uncertainties in the estimatesof hydraulic conductivity parameters in the lower layers of the model are many ordersof magnitude larger than the uncertainties of the other parameters. The uncertaintiesin these parameters dominate the model parameter space, both in the single state andjoint parameter estimates. It is also shown that prior information on these parametersstabilizes the parameter space, and leads to very small error ratios for the remainingparameters.The second stage of model calibration builds on the first stage. The hydraulicconductivity parameters in the upper model layers are the most interesting, since mostof the data is available in these layers. In the second stage, prior information is usedto specify the parameter values for the lower model layers, and the remainingparameters are estimated. The relative weighting of the hydraulic head and the 14Cdata for the joint data set is examined. The effects of possible errors in the anisotropy,source strength, and dispersivity is analyzed using the parameter space during thesecond stage of parameter estimation.2127.6 Model calibration for 10 flow parametersIn the first stage of model calibration, the model parameter space usinghydraulic head and 14C data is analyzed. Those parameters for which prior informationwill most efficiently and responsibly stabilize the model parameter set are identified.7.6.1 Parameter estimates based on the hydraulic head data setThe hydraulic head data set contains 55 hydraulic head data, 14 in theNacimiento, 18 in the Ojo Alamo, 8 in the Quaternary Alluvium, 8 in theKirtland/Fruitland, 2 in the Pictured Cliffs, and 4 in the Cliff House. The distribution ofthese data is shown in Figure 7.4. No data are available in the Lewis shale. The initialparameter estimates are listed in Table 7.4. The hydraulic conductivity parameters arethe horizontal hydraulic conductivities, in units of m/yr. For the hydraulic conductivityparameters, the initial estimates are based on the average values of the independentinformation for each of the hydrogeological units. For the head boundary parameters,the initial estimates are equal to the surface elevation at each node, in units of meters.The final parameter estimates are listed in Table 7.4, along with the standarderror of the parameter estimates and the coefficient of variation. When a different setof initial estimates were used, the final parameter estimates for the hydraulicconductivity parameters were found to be dependent on the initial estimates. The final213Table 7.4 Parameter estimates and uncertainties based on hydraulic head dataParameters Initial Estimate Final Estimate Std. Error CVKCh 1.1 1.1 21.3 498.3K 1.30e-02 1.4e-04 30.2 7.9K 7.OOe-01 8.5e-01 49.3 706.4Kk 1.40e-02 6.le-04 44.5 13.8F< 100.0 209.6 55.7 24.0K0 10.0 149.2 52.1 23.9H1 1820.0 1805.0 43.7 0.024H2 2000.0 1975.0 23.2 0.012H3 1990.0 1940.0 13.3 0.007H4 1690.0 1686.0 13.3 0.008parameter estimates for the head boundary parameters were independent of the initialestimates. Figure 7.8 shows the head distribution for the flow system using the finalparameter estimates from the head data set.For the hydraulic conductivity parameters, the standard error is given in logtransformed units. A standard error of one implies an uncertainty in the parameterestimate of one order of magnitude. The standard errors range from 21 to 55, implyingthe hydraulic conductivities are virtually non-identifiable using this data set. Theseparameter uncertainties exceed the physical bounds on hydraulic conductivity, sincethe parameter values are not constrained in the estimation procedure. For this flowmodel and boundary conditions, it could be expected that the hydraulic conductivityparameters are non-identifiable using the hydraulic head data. The coefficients of214variation (CV) is the standard error divided by the log transformed parameter value.For hydraulic conductivity parameters with an estimate near one, the log transformedestimate is small, and the CV is large. For instance, K and Kk have similar standarderrors, but the log transformed estimate for K is much smaller than Kk, and the CVof K is much larger than the CV of Kk.For the head boundary parameters, the standard errors are given in meters.The standard errors range from 13 to 43 meters, with the largest uncertainty in theestimate of the hydraulic head at the node located in the position of Chaco Canyon,node 1. From Figure 7.4, very little data exists in the vicinity of H1, and the most dataexists in the vicinity of H4, so it is reasonable that H1 has the most uncertainty and H4has the least uncertainty. The CV’s for the head parameters are the standard errordivided by the estimate, and the values of CV are relatively small. H3 has the smallestcoefficient of variation, since the standard error of H4 and H3 are nearly equal, whilethe estimate of H3 is larger than the estimate of H4. The head parameters arereasonably well defined by the data.Table 7.5 presents the fit of the simulated and the observed data, along withsome statistics describing the fit. The actual differences between the simulated andobserved data range from 2.2 m to 104.2 m. The weighted difference is the actualdifference divided by the assigned standard error for each data point. The weighteddifferences range from 0.07 to 2.69. The number of residuals greater than zero is215slightly more than the number of residuals less than zero, and the average weightedresidual is near zero. The value of the objective function at the minimum is 58.69. Thegoodness of fit can be evaluated by the s2 statistic, defined in Chapter 4 as theminimum of the objective function divided by the number of data minus the numberof parameters. If the weighted errors in the data are normally distributed and have aunit variance, the goodness of fit should be about equal to one. The goodness of fitis equal to 1.3, implying the standard errors that were assigned to the data are slightlysmaller than they should be, assuming the model is an accurate representation of theflow system.7.6.2 Parameter space based on head data setTable 7.6 lists the axes of the confidence region for the final parameterestimates using the head data set. Axis 10 is the longest axis of the confidenceregion, with a length of 690 units. The response surface is the flattest in the directionof axis 10. Axis 9 is about half as long as axis 10, and axis 8 is an order of magnitudeshorter than axis 9. All four of the longest axes are greater than 5.1 units in length,which means that the parameters which have large components from these axes havecoefficients of variation greater than 5. Examining the unit vectors, which define theorientations of the axes of the confidence region, the two longest axes only havecomponents from the hydraulic conductivity parameters. The head parameters haveno components of the unit vectors that define the orientations of the longest axes of216Table 7.5 Fit of observed and simulated head data at the parameter estimatesbased on the head dataData Identifier Observed Simulated Difference WeightedData Data Difference28n 1969.0 1936.5 32.5 0.5429n 2018.0 1952.1 65.9 1.8835n 1986.0 1958.6 27.4 0.9139n 1953.0 1898.4 54.6 1.8240n 1927.0 1898.8 28.2 1.1346n 1855.0 1824.6 30.4 0.6748n 1908.0 1864.7 43.3 0.7951n 1838.0 1835.1 2.9 0.0552n 1855.0 1869.1 -14.1 -0.2856n 1840.0 1772.6 67.4 2.69.59n 1747.0 1744.3 2.7 0.0762n 1669.0 1730.7 -61.7 -1.0372n 1660.0 1689.6 -29.6 -0.9973n 1678.0 1691.8 -13.8 -0.46190 2026.0 1949.1 76.9 1.28200 2009.0 1949.2 59.8 1,00210 1920.0 1949.1 -29.1 -0.49• 300 1932.0 1955.6 -23.6 -0.52310 1975.0 1944.1 30.9 1.54320 1912.0 1941.9 -29.9 -1.00330 1893.0 1942.7 -49.7 -2.48340 1930.0 1955.5 -25.5 -1.27360 1920.0 1943.2 -23.2 -0.42370 1930.0 1907.9 22.1 0.63380 1887.0 1910.0 -23.0 -0.66410 1935.0 1907.0 28.0 1.40430 1889.0 1924.1 -35.1 -1.75440 1876.0 1911.3 -35.3 -0.64450 1850.0 1880.7 -30.7 -0.51520 1855.0 1869.1 -14.1 -0.28610 1658.0 1729.8 -71.8 -1.31740 1654.0 1700.7 -46.7 -0.78217Table 7.5 (continued)53(1 1880.0 1873.5 6.5 0.1468ci 1663.0 1690.1 -27.1 -0.9069g 1683.0 1690.9 -7.9 -0.2670(1 1719.0 1694.4 24.6 1.23l6ci 1874.0 1868.3 5.7 0.1917p 1871.0 1864.2 6.8 0.2325(1 1956.0 1948.0 8.0 0.2726ci 1956.0 1951.2 4.8 0.168k 1862.0 1885.8 -23.8 -0.4810k 1893.0 1891.9 1.1 0.0612k 1915.0 1910.0 5.0 0.2513k 1894.0 1889.9 4.1 0.1614k 1891.0 1890.1 0.9 0.0315k 1880.0 1886.0 -6.0 -0.2020k 1935.0 1921.6 13.4 0.2227k 1787.0 1836.6 -49.6 -0.8358k 1737.0 1811.5 -74.5 -1.86ich 1898.0 1842.0 56.0 1.602ch 1852.0 1834.2 17.8 0.366pc 1938.0 1887.3 50.7 0.859ch 1731.0 1835.2 -104.2 -2.0811c 1882.0 1883.3 -1.3 -0.05l8ch 1849.0 1879.1 -30.1 -0.67Avera(1e Weighted Residual: 0.0882Number of Weighted Residuals greater than 0 : 29Number of Weighted Residuals less than 0.0 : 26Value at Objective Function Minimum = 58.69s2 = 1.3the confidence region.218Table 7.6 Axes of confidence region for parameter estimates based on head datasetAxis .1e-03 5.9e-03 i.7e-03 1.3e-02 1.9e-01 10e-01 5.le÷002.8e+O1 3.0e÷0 6.9e+02LengthsParameters U1 U2 U3 U4 U5 U6 U7 U8 U9 U10K 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 -0.87 0.50K 0.00 0.00 0.00 0.00 0.11 0.00 0.98 -0.16 0.00 0.00K 0.00 0.00 0.00 0.00 0.02 0.00 0.00 -0.02 0.50 0.87Kk 0.01 -0.01 0.02 0.00 -0.99 -0.03 0.11 -0.02 0.01 0.02K0 0.00 -0.01 -0.01 0.00 -0.02 0.71 -0.11 -0.70 -0.01 -0.01K0 0.00 0.01 0.01 0.00 0.02 -0.71 -0.11 -0.70 -0.01 -0.01H1 -0.09 0.16 -0.52 0.83 -0.02 0.00 0.00 0.00 0.00 0.00H2 -0.28 0.35 -0.71 -0.54 -0.02 0.00 0.00 0.00 0.00 0.00H3 -0.91 0.18 0.37 0.10 0.00 0.01 0.00 0.00 0.00 0.00H4 -0.30 -0.91 -0.29 -0.04 0.00 -0.01 0.00 0.00 0.00 0.00The unit vector for axis 10, the longest axis of the confidence region, showsthat the largest components come from parameter axes Kth and K, with very smallcomponents from the other parameter axes. Axis 9 is oriented in a direction betweenthe axes of KCh and but is nearly perpendicular to the axes of all other parameters.Similarly, the unit vector for axis 9, the second longest axis of the confidence region,also has the largest components from K and K, and almost no components fromthe other parameter axes. Axes 9 and 10 are perpendicular to each other, and areboth in the plane of the KCh and K axes. The plane containing axes 9 and 10contains values of the response surface very near the minimum.219The third longest axis, axis 8, has its largest components from parameter axesF< and K. The four shortest axes of the confidence region, axes 1 through 4, containmainly components from the head parameters. The confidence regions for the headparameters and the hydraulic conductivity parameters are nearly independent,because the axes which contribute to the confidence region for the head parametersare nearly perpendicular to the axes which contribute to the confidence region for thehydraulic conductivity parameters. However, the confidence region for the headparameters is not completely independent from the confidence region for the hydraulicconductivity parameters, since there are some small components of the longer sixaxes of the confidence region from the axes of the head parameters.The table of relative contributions to the CV of the parameters (Table 7.7)shows that axes 9 and 10 contribute most of the uncertainty to three of the hydraulicconductivity parameters, KCh, K,and Kk. Axis 8 contributes most of the uncertainty tothe hydraulic conductivity parameters, I< and K. The uncertainty in the headparameters comes from many different axes of the confidence region, including someof the largest axes of the confidence region. It was noted above that the headparameters have almost no projection in the direction of the six longest axes of theconfidence region. However, they do have very small projections, and these very smallprojections multiplied by the length of the long axes yield some of the uncertainty inthe head parameters.220Table 7.7 Relative contribution to the CV from each axis of the confidence ellipsoidbased on the head data setParameters U1 U2 U3 U4 U5 U6 U7 U8 U9 U10Kd, 0.00 0.00 0.00 0.00 0.00 000 0.00 0.00 0.37 0.64K1 0.00 0.00 0.00 0.00 0.00 0.00 0.53 0.43 0.02 0.03iç 0.00 0.00 0.00 0.00 0.00 0.90 0.00 0.00 0.06 0.94KR 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.06 0.94K0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.90 0.03 0.07K 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.89 0.03 0.07H1 0.00 0.00 0.03 0.27 0.03 0.00 0.04 0.08 0.48 0.08H2 0.01 0.04 0.22 0.50 0.13 0.01 0.00 0.01 0.09 0.00H3 0.38 0.03 0.17 0.05 0.00 0.25 0.00 0.04 0.04 0.03H1 0.03 0.61 0.08 0.01 0.00 0.27 0.00 0.00 0.00 0.00The parameter space has yielded the following information. First, there is alarge, multidimensional region in parameter space that has a value of the objectivefunction near the minimum. The large parameter uncertainties come from more thanone direction of parameter space. This situation is different from the multiparameterexample in Chapter 5, where large uncertainties in the parameters were present, butthese uncertainties came from only one axis of the confidence region. Second, threegroups of parameters can be identified as having nearly independent confidenceregions. The first group contains the parameters KCh, and Kk. The second groupcontains the parameters K0 and K. The third group contains the head parameters.Both the first and second group have large parameter uncertainties, but theseuncertainties are due to different axes of the confidence region.2217.6.3 Parameter estimates based on 14C data setThe model parameter space for the 14C data can be compared to the modelparameter space for the head data to determine the potential for 14C data to stabilizethe parameter estimates. The 14C data set contains 18 data, 6 in the Nacimiento and12 in the Ojo Alamo (Figure 7.6). The initial parameter estimates are chosen to beidentical to the final parameter estimates using the head data. The final parameterestimates using only the 14C data are listed in Table 7.8, along with the standard errorand coefficient of variation. Figure 7.9 shows the steady state 14C concentrationdistribution for the flow system using the final parameter estimates using only the 14Cdata. This figure shows that the 14C data exist only in the upper portion of the flowsystem, and samples from deeper formations would not result in any measurableconcentrations of 14C.For the hydraulic conductivity parameters, the standard errors for Ku,, K1, andare very large. No data are available in these units, and the 14C data is virtuallyinsensitive to these parameter values. The standard errors for I< and K are about 0.2orders of magnitude, while the standard error for Kk is about 4 orders of magnitude.For the head boundary parameters, node 1 has the largest standard error, becauseno 14C data exists in the southern end of the model domain. The remaining headnodes have standard errors between 15 and 40 meters, similar to those using headdata.222Table 7.8 Parameter estimates and uncertainties based on 14C data setParameters Initial Estimate Final Estimate Stcl. Error CVK 1.1 1.100 787.2 19088.2K1 1.44e-04 0.208 292.9 398.7K 8.52e-01 0.715 3116.4 21386.4Kk 6.lle-04 0.013 3.696 1.893iç 209.6 10.770 0.204 0.203K 149.2 0.398 0.167 0.363H1 1805.0 1983.0 4153.1 2.218H2 1975.0 1991.0 33.2 0.017H3 1940.0 1997.0 15.0 0.008H4 1686.0 1714.0 39.8 0.023Table 7.9 shows the fit of the simulated and the observed 14C data. The actualdifferences between the simulated and observed data 14C concentrations range from0.05 to 11.1 pmc. The weighted differences ranged from 0.10 to 3.08. The residualsgreater than zero outnumber residuals less than zero, and the average weightedresidual is greater than zero. The value of the objective function at the minimum is40.5. The goodness of fit is 5.0, larger than the goodness of fit for the head data.Again, the assigned standard errors are smaller than they should be, according to howwell the data fit the model.223Table 7.9 Fit of observed and simulated 14C data using parameter estimates basedon 14C data setData Identifier Observed Data Simulated Data Difference Weighted Difference28n-c 51.20 52.59 -1.39 -0.1539n-c 12.60 11.37 1.23 0.5440n-c 19.90 17.63 2.27 0.9148n-c 9.20 19.07 -9.87 -2.24Bin-c 1.40 0.03 1.37 2.7482n-c 4.50 3.94 0.56 0.14200-c 37.90 32.97 4.93 0.63220-c 10.16 8.71 1.45 0.29230-c 50.50 54.40 -3.90 -0.52320-c 27.20 38.61 -11.41 -3.08330-c 40.30 29.78 10.52 2.63340-c 5.52 5.57 -0.05 -0.10380-c 8.33 11 .25 -2.92 -1.27520-c 8.14 5.79 2.35 0.62570-c 2.79 0.48 2.31 0.62580-c 3.96 2.89 1.07 0.45830-c 8.85 4.28 4.57 2.54840-c 4.70 5.77 -1.07 -0.76Average Weighted Residual 0.096Number of Weighted Residuals greater than 0: 11Number of Weighted Residuals less than 0.0: 7Value at Objective Function Minimum = 40.49s2 = 5.02247.6.4 Parameter space based on 14C data setTable 7.10 lists the axes of the confidence region for the final parameterestimates using the ‘4C data set. The longest axis is axis 10; the second longest isaxis 9. Both of these axes have significant components from the parameter axes KChand K, a small component from K1, and almost no components from the otherparameters. Axis 8 is 18 units long, and is aligned closely with the parameter axis K1.Axis 7 is 1.6 units long, and is aligned most closely with the axis of Kk. All other axesTable 7.10 Axes of confidence region for parameter estimates based on 14C datasetAxis 9.2e-04 8.Oe-03 1.le-02 5.2e-02 1.3e-01 2.9e-O1 1.6e÷O0 1.8e+01 4.7e+02 2.9e+04LengthsParameters U1 U2 U3 U4 U5 U6 U7 U8 U9 U10Kth 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.03 0.75 0.67K1 0.00 0.00 0.00 -0.04 0.04 0.00 0.05 -1.00 0.06 -0.01K 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -0.05 -0.66 0.75K 0.03 -0.01 0.01 -0.25 0.22 -0.17 -0.93 -0.03 0.00 0.00K0 -0.05 -0.02 0.01 0.65 0.69 0.31 -0.07 0.00 0.00 0.00K0 -0.03 0.04 -0.04 0.41 0.03 -0.91 0.07 -0.01 0.00 0.00H1 0.07 -0.03 0.02 -0.59 0.69 -0.22 0.35 0.06 0.00 0.00H2 0.45 -0.66 -0.60 0.04 -0.02 0.01 0.00 0.00 0.00 0.00H3 0.89 0.30 0.35 0.08 -0.01 0.00 0.00 0.00 0.00 0.00H4 0.05 0.69 -0.72 -0.03 0.04 0.05 -0.01 0.00 0.00 0.00225are less than one unit long. The table of relative contributions to the CV of theparameters (Table 7.11) shows that axis 10 contributes most of the uncertainty tothree of the hydraulic conductivity parameters, Kd,, K0,and K1, and the parameter H1.Table 7.11 Relative contribution to the CV from each axis of the confidence regionbased on 14C data setParameters U1 U2 U3 U4 U5 U6 U7 U8 U9 U10Kct, 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00K1 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.99K0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00Kk 0.00 0.00 0.00 0.00 0.00 0.00 0.63 0.07 0.29 0.02K0 0.00 0.00 0.00 0.03 0.21 0.21 0.31 0.04 0.19 0.02K 0.00 0.00 0.00 0.00 0.00 0.56 0.08 0.34 0.01 0.00H1 0.00 0.00 0.00 0.00 0.00 0.00 0.07 0.25 0.11 0.58H2 0.00 0.10 0.16 0.02 0.02 0.01 0.11 0.52 0.06 0.00H3 0.01 0.10 0.26 0.33 0.04 0.01 0.15 0.03 0.00 0.07H1 0.00 0.06 0.12 0.00 0.05 0.42 0.12 0.02 0.00 0.22The confidence region using 14C data can be compared to the confidenceregion using the head data by calculating angles of interaction. For the longest axes(axis 10 for ‘4C data; axis 10 for head data), the angle of interaction is 8.9°, indicatingthat the orientations of the axes are not much different. For the second longest axes(axis 9 for both data sets), the angle of interaction is also 8.9°. Recall that the longestaxes contain mainly uncertainties from the hydraulic conductivity parameters from thelower layers of the model. Because the axes are very long, even this small angle of226interaction reduces the uncertainty of the parameters using the joint data set.However, the uncertainty of the parameters is still be large, and the maximumuncertainty still occurs along this direction of parameter space. It is more difficult topredict the reduction in uncertainty for the other model parameters.7.6.5 Joint parameter estimationJoint parameter estimation using both hydraulic head and 14C data is used toestimate the model parameters. For joint estimation, the weights for each data set arerequired. The weights can be chosen using any of the methods presented in Chapter6. Table 7.12 contains the parameter estimates and uncertainties for all four weightingcriteria. All four criteria weight the head data more heavily than the concentration data.The head weights all have Wh = 1 .0, and the concentration weights range from w0 =0.1 to w = 0.52.Table 7.12 also contains some additional statistics about each set of parameterestimates, including the average CV, the CN, the volume of the confidence region,and the length of the longest axis of the confidence region. It also contains statisticsabout the fit of the observed and simulated data at each set of parameter estimates.Sh is the value of the head portion of the objective function evaluated at the set ofparameter estimates, and S is the value of the concentration portion of the objectivefunction evaluated at the same estimates. The h2 statistic is Sh/nh, and the s2 statistic227r\)CD)-‘CD)CDCD0’ CD Cr,D a C D CD CD ::i.D CD 0 —‘ 0 D -I a 1) C,) CD C C,) D CD -.‘0 CH 0 0 CDHeadDataOntyConcentrationDataOntyointDataSet;MINCNointDataSet;RESIDtointDataSet;MINLENJointDataSet;MINVOLWh1.00.01.01.01.01.0w0.01.00.100.190.240.54ParametersEstimateStdCoeff.EstimateStdCoeff.EstimateStdCoeft.Var.EstimateStdCoeff.EstimateStdCoeff.EstimateStdCoeff.Var.Var.Var.Var.VarK1.lOe÷0021.29498.31.100787.219088.21.10e+004.49108.7001.lOe+004.50108.9501.lOe+004.51109.5001.100+004.29104.200K1,44e-0430.197.90.208292.9396.71.220-034.371.5001.39e-034.731.6571.52e-034.671.6552.41e-034.221.613K.,8.52e-0149.33706.40.7153116.421386.46.99e-015.8037.2407.OOe-015.2134.7507.OOe-014.7430.5707.OOe-018.9457.830K6.lle-0444.4613.80.0133.6961.8931.64e-030.300.1061.57o-030.260.0931.55e-030.260.0931.47e-030.350.123K2.lOe+0255.7324.010.7700.2040.2035.Ole+000.220.3085.67e÷000.180.2346.02e÷000.150.1987.32e+000.110.122K1.49e÷0252.0523.90.3980.1670.3631.64e-010.180.2241.78e-010.140.1881.88e-010.130.1752.16e-010.100.151H,1805.043.720.0241983.04153.12.2181803.028.420.0161804.028.880.0161804.029.740.0161800.033.910.019H1975.023.200.0121991.033.20.0171971.016.140.0081972.014.640.0071973.013.810.0071982.012.050.006H1940.013.330.0071997.015.00.0081952.08.270.0041956.07.770.0041959.07.450.0041969.06.780.003H1686.013.250.0081714.039.80.0231679.010.690.0061678.010.800.0061677.011.100.0071676.012.000.007AverageCV127.444.le÷0314.81114.59114.22316.407ConditionNumber1.7e+052.8e+073.9e+044.3e+044.5e÷045.Ie+04VolumeofRegion4.73e-O5.37e-051.25e-098.29e-104.86e-12.51e-10LengthofMaxAxis6902.90e+04131.2121.3111.9117.8S,58.70.062.763.866.070.7S,0.040.5100.797.585.272.1s1.060.001.141.161.201.2820.002.225.595.414.734.00is SJnc, where nh and n,, are the number of head and 14C data respectively. The s2statistics give a measure of the fit between the observed and simulated data, and areused to compare fits among the different weighting methods. For the purpose ofcomparing fits, the number of parameters being estimated is not contained in the s2statistic.It is important to note that the fit for each individual data set using the jointestimates is somewhat worse than the fit for each data set using the estimates fromthat data set. Each data set leads to a different set of parameter estimates, and thejoint parameter estimates reflect the trade-off between the two data sets. For instance,using the parameter estimates based on the head data only, Sh is 58.7, while for thejoint parameter estimates using the MINLEN criterion, Sh is 66.0. Similarly, using theparameter estimates based on the 14C data only, S,, is 40.5, while for the jointparameter estimates using the MINLEN criterion, S,, is 85.2. The difference betweenSh using the head data set and Sh using the joint data set is smaller than thedifference between S,, using the 14C data set and S using the joint data set, reflectingthe fact that the head data is weighted more heavily than the 14C data. Examining thes2 statistics, the ‘4C data always fits worse than the head data, for all sets of weights.However, as the relative 14C weights increase, the fit of the 14C data becomes betterand the fit of the head data becomes worse.The parameter estimates and uncertainties for all four weighting criteria are229similar, though definite patterns emerge. The standard errors and CV’s for allparameters are somewhat reduced from those using only the head data. Theparameters KCh, K1 and K. all have standard errors near 5 orders of magnitude.Though these standard errors are reduced from those produced using head ordata alone, they are still significant. The parameters Kk, K0 and K all have standarderrors approximately between 0.1 and 0.3, much reduced from those using only headdata. The uncertainty in Kk is significantly reduced from that using either data setalone, though the uncertainty in Kk increases as relative weighting of the 14C dataincreases. For K and K0, the uncertainties are nearly the same or somewhat largerthan those using the 14C data alone. Note that the uncertainties in K and I< decreaseas the relative weighting of the concentration data increases. The head parametershave standard errors between 7 and 33 meters, reduced by nearly a factor of twofrom those using either data set individually.Comparing the different weighting criteria, the MINCN, RESID, and MINLENcriteria all have relatively small concentration weights. The parameter estimates areslightly different for each set of weights. The smallest average CV is produced by theMINLEN criterion, though the average CV’s for the RESID criterion and MINCNcriterion are not much larger. Both the MINLEN and the RESID criterion seem to strikea good balance between minimizing the condition number and minimizing the volumeof the confidence region. It must be noted that the average CV’S are dominated by thelarge CV’S for parameters K and230Figure 7.10 shows the head distribution for the flow system using the jointparameter estimates using the MINLEN criterion. Comparing Figure 7.10 with Figure7.8, the head distribution using only head data, the joint data set seems to produceparameter estimates which allow more flow through the lower layers of the system.The joint data set results in larger conductivities for the aquitards, which allows greaterleakage into the lower model layers. Figure 7.11 shows the concentration distributionfor the flow system using the same estimates.The four weighting criteria all weight the head data more than the concentrationdata. For the parameters with the largest uncertainties, these weights reduce theparameter uncertainties as much as possible. However, for the parameters F< and K,the 14C data results in much better parameter estimates than the head data. Weightingthe concentration data less than the head data does not result in the smallest possibleuncertainties in these parameters. As the relative weight of the concentration datadecreases, the fit of the 14C data becomes worse. The estimates of those parameterswhich the 14C data can estimate well (K0 and K in this case) becomes more uncertain.7.6.6 Parameter space based on joint data setTable 7.13 lists the axes of the confidence region for the final parameterestimates using the joint data set for the MINLEN criterion. The confidence region hasmany similarities to the confidence region for the head data set (Table 7.6). Axes 9231and 10 are the longest axes of the confidence region. Both have significantcomponents from the parameter axes K and K, but small components from theother parameter axes. The third longest axis, axis 8, is nearly parallel to parameteraxes K1. The fourth and fifth longest axes, axes 7 and 6, have their largestcomponents from parameter axes K0 and K.Table 7.13 Axes of confidence region for parameter estimates based on joint dataset at MINLEN weightsAxis 2.6e-03 4.3e-03 6.2e-03 1.2e-02 4.4e-02 1.6e-01 2.le-01 2.9e-01 1.le+02 1.3e+01LengthsParameters U1 U2 U3 U4 U5 U6 U7 U8 U9 U10Kd, 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -0.97 0.25K1 0.00 -0.01 0.00 -0.04 0.14 -0.01 0.08 -0.99 0.01 0.06K0 0.00 0.00 0.00 0.00 0.01 0.00 0.01 -0.06 -0.25 -0.97Kk -0.01 -0.10 0.01 -0.11 0.97 0.10 0.05 0.14 0.00 0.00K0 -0.01 0.00 0.00 0.00 -0.10 0.47 0.87 0.05 0.00 0.00K 0.01 -0.01 -0.01 0.00 -0.05 0.88 -0.48 -0.05 0.00 0.00H1 0.01 0.02 -0.03 -0.99 -0.12 -0.01 -0.01 0.02 0.00 0.00H2 0.18 0.97 -0.13 0.01 0.10 0.01 0.01 0.00 0.00 0.00H3 0.95 -0.14 0.27 0.00 -0.01 -0.01 0.01 0.00 0.00 0.00H4 0.25 -0.18 -0.95 0.02 0.00 -0.01 0.01 0.00 0.00 0.00The table of relative contributions to the CV of the parameters (Table 7.14)shows that axis 9 contributes most of the uncertainty to four of the hydraulicconductivity parameters, KCh, K1, K,and Kk. Axes 6 and 7 contribute most of the232uncertainty to the other two other hydraulic conductivity parameters, K0 and K. Theuncertainty in the head parameters comes from many different axes of the confidenceregion.Table 7.14 Relative contribution to the CV from each axis of the confidence ellipsoidbased on the joint data set at M[NLEN weightsParameters U1 U2 U3 U4 U5 U6 U7 U8 U9 U10KCh 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00K1 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.03 0.79 0.19K0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.84 0.16Kk 0.00 0.00 0.00 0.00 0.22 0.03 0.02 0.20 0.36 0.191<0 0.00 0.00 0.00 0.00 0.00 0.15 0.84 0.01 0.00 0.00K0 0.00 0.00 0.00 0.00 0.00 0.67 0.32 0.01 0.00 0.00H1 0.00 0.00 0.00 0.49 0.10 0.00 0.01 0.15 0.16 0.09H2 0.01 0.37 0.01 0.00 0.43 0.11 0.06 0.01 0.00 0.00H3 0.43 0.03 0.21 0.00 0.03 0.06 0.24 0.01 0.00 0.01H4 0.01 0.01 0.84 0.00 0.00 0.08 0.05 0.00 0.01 0.00The error ratio matrix (Table 7.15) shows smallest error ratios result from priorinformation on KCh, while the second smallest error ratios result from prior informationon K. Prior information on any of the parameters KCh, K1, or Kk result in verysmall error ratios for the remaining parameters. Errors in any other parameters leadto very large error ratios for the parameters KCh, and K1. Since these parametershave large uncertainties, the errors due to the error ratio may be within the parameteruncertainties. Regardless, prior information on the parameters KCh, K0, and K1 are the233r\) ()Om—-s(DO_5_s00 3 -5 x -.‘ 0 —5 0 D 3 CD -‘ CD -5 CD 0 3 CD (I) (1) CD 0 0 D z m z- cr CD C,,ParameterswithpriorinformationEstimatedKc,K1KKkK0K0H1H2H3H4ParametersKCh1.00-58.603.28-700.286.4410.022626.75-281.21-683.161437.06K1-0.011.00-0.049.38-0.040.00-35.203.779.15-19.26K0.26-14.971.00-178.886.03.80670.98-71.83-174.51367.08Kk0.000.030.001.000.000.00-1.330.140.35-0.73K00.000.000.00-0.011.000.000.030.00-0.010.01K0.000.000.000.000.001.000.000.000.000.00H,0.000.000.00-0.040.000.001.00-0.02-0.040.09H20.000.000.000.000.000.000.001.000.000.00H30.000.000.000.000.000.000.000.001.000.00H40.000.000.000.000.000.000.010.000.001.00most responsible since they lead to the smallest error ratios.Based on examining the parameter space, three groups of parameters can bedefined. The first group contains parameter Kth, K, K1 and Kk. For the joint data set,prior information on Kth is the most efficient in reducing the overall parameteruncertainty, and leads to the smallest linearized error ratios. Prior information on anyof the parameters in this group will reduce the uncertainty for all the parameters in thegroup. The second group contains the parameters 1<0 and K. Prior information onthese parameters will not reduce the uncertainty of any of the parameters in the firstgroup. Prior information on K0 or K will mainly reduce the uncertainty in these twoparameters. The third group contains the head parameters. Prior information on thehead parameters will not reduce the uncertainty in any of the hydraulic conductivityparameters significantly. The linearized error ratios show that the head parameters arethe least responsible parameters for prior information.7.6.7 Discussion of results for 10 parameter systemEstimating all 10 flow parameters in the model resulted in very largeuncertainties for the hydraulic conductivities in the lower model layers. These largeuncertainties dominated the model parameter space, and thus dominated the way theparameter space criteria selected the weights for joint parameter estimation. Priorinformation on these parameters with large uncertainties will efficiently stabilize the235parameter set. From the error ratio matrix, errors in the prior information for any ofthese parameters will not result in significant errors in the parameter estimates of theother model parameters.Because the large uncertainties in the hydraulic conductivities of the lowermodel layers dominated model parameter space, the weighting criteria focused onreducing these large uncertainties. However, if the modeller were interested inestimating the hydraulic conductivity of the upper model layers as well as possible, theresults would be disappointing, since the weighting criteria did not lead to the bestestimates of these parameters. Since the uncertainties in the hydraulic conductivitiesof the upper model layers was small compared to the lower model layers, reducingthe uncertainties in the upper model layers was not a priority for the weighting criteria.In fact, the weighting criteria that weighted the head data more than the 14C dataincreased the uncertainties in the upper model layers.For the second stage of parameter estimation, we wish to focus on estimatingthe hydraulic conductivities of the upper model layers as well as possible, given thedata available. Based on the first stage, it is efficient and responsible to use priorinformation to define the parameter values for the lower four hydraulic conductivityparameters, KCh, K, and Kk. These four parameters are specified at the parametervalues from joint parameter estimation. Specifying the parameter values of the lowermodel layers will not significantly affect the parameter estimates for the upper model236layers.7.7 Model calibration for 6 flow parametersFor the second stage of model calibration, only six model parameter areestimated. These parameters are the two hydraulic conductivity parameters, iç andK, and the four head parameters. The same model and boundary conditions areused. The other four hydraulic conductivity parameters are specified at the valuesfrom the joint parameter estimates, which are very close to the prior parameterestimates. The model parameter space using hydraulic head and 14C data is analyzed.The weighting criteria for the joint parameter estimates is the main focus. Afterestimation of the parameters using the joint data set is disussed, the effect of errorsin the other model parameters is examined.7.7.1 Parameter estimates and confidence region based on hydraulic head datasetThe parameter estimates for this reduced parameter set using only hydraulichead data are nearly identical to the parameter estimates for the full parameter setusing hydraulic head data. Table 7.16 contains the parameter estimates, standarderrors, and CV’S for the 6 parameters using only hydraulic head data. The standarderrors and CV’s are somewhat reduced from the full parameter set. With the full237parameter set, some of the longest axes of the parameter confidence region containedpart of the uncertainty in these parameters. These long axes are now eliminated, andtherefore the uncertainties in these parameters is somewhat reduced. However, thestandard errors for the hydraulic conductivity parameters are still extremely large,suggesting these parameters are virtually non-identifiable using head data. Thestandard errors for the head parameters range from 9.8 to 25 meters.Table 7.17 lists the axes of the confidence region for the 6 parameter set usingthe head data. Axes 5 and 6 are the longest axes, and have components only fromthe parameter axes K0 and K. Axes 1 through 4 contain components only from thehead parameters. Again, the confidence regions for the hydraulic conductivityparameters and head parameters are nearly independent of each other. Note thataxes 1 through 4 of the 6 parameter system are very similar to axes 1 through 4 ofthe 10 parameter system (Table 7.6). Also, axes 5 and 6 of the 6 parameter systemare similar to axes 6 and 8 of the 10 parameter system. These similarities show thatthe confidence region for the parameters in the 6 parameter system is nearly thesame as that for the 10 parameter system. The table of relative contributions to theCV of the parameters (Table 7.18) shows that axis 6 contributes most of theuncertainty to the two hydraulic conductivity parameters. The four smallest axescontribute most of the uncertainty to the head parameters.238CA)(0(Do)3CD-I.CQCDC)‘D.CD Cl)1) .D(DOCnCn-. DCD ?C/) -4,0 0 D a C)) CD C/) CD a. 0 - 0 cCD 0)HeadDataOnlyConcentrationDataOnlyJointDataSet;MINCNJointDataSet;RESIDJointDataSet;MINLENJointDataSet;MINVOLWh1.00.01.01.00.600.19w,0.01.00.170.321.01.0ParametersEstimateStdCoeff,EstimateStdCoeff,EstimateStdCoeff,EstimateStdCoeff.EstimateStdCoeft.EstimateStdCoeff.Var.Var.Var.Var.Var.VarK209.33.7914,56010.7500.1370.1355.3680.1480.2025,3610.1020.1407.6680.0730.0829.3360.0630.065Kn149.33.7715.5300.3970.1300.2760.2160.1190.1780.2340.0890.1410.2580,0660.1120.3480.0620.135H,1805.025.040.0141836.0352.080.1921817.020.010.0111812.021.280.0121787.032.360.0181766.048.370.027H11975.019.520.0101991.016.950.0091968.011.440.0061975.09.440.0051995.07.900.0041994.06.500.003H,1940.09.830.0051997.09.480.0051960.06.900.0041969.05.960.0031987.04.680.0021988.04.180.002H,1686.013.350.0081714.033.180.0191678.010.380.0061677.011.090.0071674.012.520.0071699.013.510.008AverageCV5.0210.1060.0680.0510.0380.040ConditionNumber4716.0321.089.592.998.6135.9VolumeofRegion1.55e-01.Ole-13.02e-111.lOe-112.63e-1:1.96e-12LengthofMaxAxis19.10.2510.2190.1690.1260.143Sn58.70.06.4.473.589.2101.8S0.040,594.173.262.556.7sr,’1.060.001.171.331.621.85s320.002.225,224.073.473.15Table 7.17 Axes of confidence region for parameter estimates based on head data;6 parameter systemAxis Lengths 4.le-03 6.le-03 7.le-03 1.4e-02 3.4e+00 1.9e+01Parameters U1 U2 U3 U4 U5 U6K0 0.00 0.00 0.00 0.00 -0.73 -0.68K0 0.00 0.00 0.00 0.00 0.68 -0.73H1 -0.06 -0.08 -0.49 -0.87 0.00 0.00H2 -0.23 -0.23 -0.81 0.49 0.00 0.00H3 -0.94 -0.19 0.27 -0.08 0.00 0.00H4 -0.25 0.95 -0.19 0.03 0.00 0.00Table 7.18 Relative contribution to the CV from each axis of the confidence regionbased on head data; 6 parameter systemParameters U1 U2 U3 U4 U5 U6K0 0.00 0.00 0.00 0.00 0.03 0.97K0 0.00 0.00 0.00 0.00 0.03 0.97H1 0.00 0.00 0.07 0.88 0.00 0.04H2 0.01 0.03 0.40 0.55 0.00 0.01.H3 0.68 0.06 0.17 0.05 0.04 0.00H4 0.02 0.65 0.03 0.00 0.28 0.027.7.2 Parameter estimates and confidence region based on 14C data setTable 7.16 also contains the parameter estimates and uncertainties for the 14Cdata set when estimating only 6 parameters. The parameter estimates for this reduced240parameter set using only 14C data are nearly identical to the parameter estimates forthe full parameter set using the same data. The standard errors and CV’S are againsomewhat reduced from the full parameter set, for the same reasons explained above.The standard errors for the hydraulic conductivity parameters are near 0.13 orders ofmagnitude. The largest standard error is for parameter H1, due to the scarcity of datain the southern end of the model.Table 7.19 lists the axes of the confidence region forthe 6 parameter set usingthe 14C data. Axis 6 is the longest axis. Axes 4, 5 and 6 have the largest componentsfrom parameter axes K0, K and H1. Axes 1 through 4 have the largest componentsfrom H1, H2, and H3. The six axes of the 6 parameter system follow the same patternas axes 1 through 6 of the 10 parameter system (Table 7.10). The table of relativecontributions to the CV of the parameters (Table 7.20) shows that axes 5 and 6contributes most of the uncertainty to K0, K, H1, and H4, while axes two and threecontribute most of the uncertainty to H2 and H3.The confidence region using 14C data can be compared to the confidenceregion using the head data by calculating the angle of interaction. For the longest axes(axis 6 for 14C data; axis 6 for head data), the angle of interaction is 37°. The axeshave quite different orientations. However, using the head data, the longest axis is 19units long, while for 14C data, the longest axis is 0.25 units long. Though the axeshave different orientations, the axis for 14C data is much shorter than the axis for head241Table 7.19 Axes of confidence region based on 14C data; 6 parameter systemAxis Lengths 8.le-04 6.5e-03 9.Oe-03 7.Oe-02 1.3e-01 2.5e-01Parameters U1 U2 U3 U4 U5 U6K0 -0.01 0.03 -0.04 0.80 -0.54 -0.25K0 0.00 -0.02 0.04 0.44 0.25 0.86H1 0.02 0.01 0.01 -0.40 -0.80 0.44H2 0.43 0.78 0.46 0.01 0.01 -0.01H3 0.89 -0.29 -0.34 0.01 0.01 0.00H4 0.14 -0.55 0.82 0.02 -0.03 -0.05Table 7.20 Relative contribution to the CV from each axis of the confidence ellipsoidbased on the 14C data; 6 parameter systemParameters U1 U2 U3 U4 U5 U6K0 0.00 0.00 0.00 0.26 0.41 0.33K0 0.00 0.00 0.00 0.02 0.02 0.96H1 0.00 0.00 0.00 0.03 0.44 0.53H2 0.00 0.50 0.34 0.02 0.05 0.09H3 0.03 0.22 0.63 0.04 0.07 0.01H4 0.00 0.05 0.21 0.01 0.08 0.65data. The parameters are much better defined using the 14C data; thus the joint dataset may not improve the parameter estimates significantly over those using 14C dataalone. Remember that the longest axes contribute mainly to the uncertainty in thehydraulic conductivity parameters.2427.7.3 Joint parameter estimation for 6 parameter systemJoint parameter estimation using both hydraulic head and 14C data is used toestimate the 6 model parameters. For joint estimation the weights for each data setcan be chosen using any of the methods presented in Chapter 6. Before presentingthe weights calculated using each method, the issue of weighting for this parameterset is discussed in a conceptual manner.Table 7.16 contains the parameter estimates, uncertainties, and other statisticsabout the parameter estimates, including information about the fit for each data set.The head data produce better fits than the 14C data at all sets of weights. However,even with the worse fit, the 14C data estimate the two hydraulic conductivityparameters much better than the head data, and the 14C data estimate the H2 and H3about as well as the head data set. When the head data is weighted more heavilythan the 14C data, the fit of the calculated head data to the observed head data isbetter, since the parameter estimates are closer to those using only head data. At thesame time, the fit for the 14C data is worse, since the parameter estimates are fartherfrom the estimates using only 14C data. Since the calculated uncertainty depends onthe fit of the data (see equations (4.2) and (4.3)), those parameters which the 14C dataestimate well (K0 and K in this case) is poorly served by weighting the head datamore than the 14C data. Conversely, if the 14C data are weighted more heavily thanthe head data, the fit cf the calculated 14C data to the observed 14C data is better,243since the parameter estimates are closer to those using only 14C data. Thoseparameters which the 14C data estimate well, K0 and K, benefit by weighting the 14Cdata more than the head data. It makes sense to weight the 14C data more heavilythan the head data in this case, since the ‘4C data produces better estimates of thehydraulic conductivity parameters than the head data.Table 7.16 contains the parameter estimates and uncertainties for all fourweighting criteria. The MINCN and RESID criteria weight the head data more heavilythan the 140 data. The RESID criterion produces weights which are somewhat differentfrom the weights produced from the 10 parameter set. The RESID criterion weightsthe head data more heavily because Sh2 is always smaller than s2. The MINLEN andMINVOL criteria weight the 14C data more heavily than the head data. As discussedabove, this weighting makes sense conceptually, since the 14C data produce betterestimates than the head data for the most uncertain parameters.The parameter estimates and parameter space change as the weights change.For K and K0, the parameter estimates become larger as the 14C weight increases.The uncertainties in these parameters also become decrease as the 14C weightincreases. For parameters H1 and H4, the parameter uncertainty increases as the 14Cweight increases. These parameters are more uncertain using 14C data than headdata, so this increase in parameter uncertainty is not surprising. For parameters H2and H3, the parameter uncertainties decrease slightly as the 14C weight increases.244Though individual parameter uncertainties move different directions as the 14Cweight increases, the overall parameter uncertainty generally decreases as the ‘4Cweight increases. The minimum overall parameter uncertainty, as measured by theaverage CV, is observed at the parameter estimates using the MINLEN weightingcriterion. Weighting the head data more than the ‘4C data, as required by the MINCVand RESID criteria, does not result in the lowest uncertainty joint parameter estimates.7.7.4 Parameter space based on joint data setIf more reliable parameter estimates are required, the parameter space for the6 parameter set can be examined to determine the most efficient and responsibleparameters for prior information. Table 7.21 lists the axes of the confidence region forthe final parameter estimates for the joint data set using the MINLEN parameterestimates. The longest axis, axis 6, has the largest component from parameter K,and the second largest component from axis K0. Prior information on these twoparameters would be the most efficient in reducing the overall parameter uncertainty.Table 7.22 is the error ratio matrix for this parameter set. Prior information onK results in the smallest error ratios. Prior information on K0 also results in relativelysmall error ratios, as does prior information on H1. Prior information on the other headparameters lead to larger error ratios for the hydraulic conductivity parameters. K isthe most responsible parameter for prior information, while K0 and H1 are also good245Table 7.21 Axes of confidence region based on joint data set at MINLEN weights;6 parameter systemAxis Lengths 1 .29e-03 3.64e-03 4.47e-03 1 .76e-02 1 .26e-01 5.50e-02Parameters U1 U2 U3 U4 U5 U6I< -0.03 0.02 0.00 0.05 -0.53 -0.85K -0.01 0.02 0.05 0.05 0.85 -0.53H1 0.04 -0.03 -0.02 -1.00 0.02 -0.07H2 0.24 -0.95 0.22 0.04 -0.01 -0.03H3 0.96 0.20 -0.19 0.04 -0.01 -0.03H4 0.14 0.25 0.96 -0.02 -0.05 0.02Table 7.22 Error ratio matrix based on joint data set at MINLEN weights; 6parameter systemParameters with prior informationEstimated K0 K0 H1 H2 H3 H4ParametersI< 1.00 -0.58 -0.51 6.71 17.91 5.07K0 -1.08 1.00 0.81 -1.72 -4.29 9.24H1 -0.02 0.02 1.00 -0.46 -0.32 -0.26H2 0.01 -0.01 0.00 1.00 0.08 0.06H3 0.01 -0.01 -0.01 0.03 1.00 0.07H4 0.06 -0.05 -0.04 0.23 0.75 1.00choices for prior information. Prior information on the other head parameters wouldnot be as responsible if the modeller were interested in the estimates of I< and K.2467.7.5 Effect of errors in other parametersIn calibrating the above model, several parameters were not estimated. Theseparameters are the effective porosity in each unit, the anisotropy for each unit, theinitial concentration for 14C at the boundary, and the longitudinal and transversedispersivities. These parameters were assigned a specified value. It is important todetermine whether errors in the specified values of these parameters influence theestimated parameters significantly. For this purpose, the linearized error ratio matrixcan be computed for these transport parameters.The resulting error ratio matrix is shown in Table 7.23. The parameters P0 andP, are the effective porosity of the Ojo Alamo and Nacimiento respectively, and theparameters An0 and An are the anisotropy of the Ojo Alamo and Nacimiento.Theparameters that have significant error ratios are the hydraulic conductivity parametersas a result of errors in c, c, and Ann. Errors in the boundary 14C concentration maybe magnified over three times when estimating the hydraulic conductivity parameters.Errors in cL lead to largest error ratios for K0 and K, while errors in xT lead to verysmall error ratios for K than I<. The estimates of the hydraulic conductivityparameters are very sensitive to errors in the estimate of c. This sensitivity isproblematic, since prior information on aL is difficult to obtain. The parameterestimates are insensitive to errors in the anisotropy of K0, but more sensitive to theanisotropy of K. These error ratios are consistent, since in the Ojo Alamo, the flow247Table 7.23 Error ratio matrix for transport parameters using parameter space basedon joint data set at MINLEN weights; 6 parameter setTransport ParametersEstimated P0 P, c An0 AnParametersK0 0.55 0.52 -3.6 -9.0 -0.30 -0.01 0.61K -0.79 -0.71 3.3 11 .5 0.86 0.01 -1.37H1 0.00 0.00 -0.01 0.02 0.00 0.00 0.00H2 0.00 0.00 0.01 -0.01 0.00 0.00 0.00H3 0.00 0.00 0.00 0.00 0.00 0.00 0.01H4 0.00 0.00 -0.02 0.00 0.01 0.00 0.01direction is dominantly along the dip of the unit, and therefore not very sensitive to theanisotropy of K0. In the Nacimiento, there is a small component of vertical flow, andthe anisotropy of the unit does affect the flow field. Based on these linearized errorratios, estimating both and XL along with the flow parameters for this model wouldbe more responsible than simply using prior information on these transportparameters. Using prior information for the other transport parameters is responsible,since the linearized error ratios are small.Based on the above error ratio analysis, parameter estimation for the flowsystem is performed for 8 parameters: the six flow parameters and and Forthe final parameter estimate is 79.5 pmc. For aL, the final parameter estimates is 350m. Both of these parameter estimates are close to their prior values, so the estimatesfor the flow parameters are similar to the parameter estimates using prior information248on c1 and c. The linearized error ratios show the potential for large errors whenusing prior information on cjr and aL, but these prior estimates are close to theestimates using the calibration data, so no large errors occur.7.8 Summary of model calibration for the San Juan BasinA cross section model for the San Juan basin has been developed. Hydraulichead data, 14C data, and independent information on parameter values have beenused to calibrate this model in an efficient and responsible manner. The model hasbeen calibrated in two stages.In the first stage, 10 flow parameters were estimated, the hydraulic conductivityof the 6 model layers and the head boundary values at four nodal points. The jointestimates did reduce the uncertainty in the model parameters significantly over thehead data set. However, estimating all 10 parameters resulted in large uncertaintiesfor the hydraulic conductivities in the lower model layers in both single state and jointparameter estimates. These large uncertainties dominated the model parameterspace, and thus dominated the way the parameter space criteria selected the weightsfor joint parameter estimation. The criteria weighted the head data more heavily thanthe concentration data. These weights reduced the uncertainty in the overallparameter estimates, but did not produce the best estimates for all model parameters.Prior information on the parameters in the lower model layers was found to be most249efficient in stabilizing the overall parameter set. Errors in the prior information for anyof these parameters did not result in significant errors in the parameter estimates ofthe other model parameters.In the second stage of model calibration, the focus was on obtaining the bestpossible estimates for the hydraulic conductivities in the upper model layers. Priorinformation was used to define the parameter values for the lower four hydraulicconductivity parameters, and only six model parameters were estimated: two hydraulicconductivity parameters and the four head parameters. The concentration dataresulted in much better estimates of the hydraulic conductivity parameters than thehead data set. For this parameter set, the MINLEN and MIN VOL parameter spacemethods weighted the concentration data more heavily than the head data. Theseweights made sense conceptually, and resulted in the best estimates for the hydraulicconductivities in the upper model layers, and the minimum overall parameteruncertainties. The RESID and MINCN criteria weighted the head data more heavilythan the 14C data, a weighting scheme which did not result in the best parameterestimates. The error ratios for the parameters which were specified were generallysmall, with the exception of the initial concentration of 14C and the longitudingaldispe rsivity.250Figure 7.1 Plan view of the San Juan Basin showing outcrops of major units and thelocation of the modelled cross section. The white in the center of the basinrepresents the outcrop of the Nacimiento Formation. The shaded areasrepresent the outcrop area of the Ojo Alamo SS, Kirtland/Fruitland Shale,and Pictured Cliffs SS.HSOUTH4RN._L__Z____-LETf[kNt________—1—lgtonf•1TIQt(*t0L.1)E7r\_ -RI( ARRIB1\ )_.JS./ICKrnMY]251OutcropofOjoAlamoSS-I, C -‘ CD t\)3Grn)CDQD(DCD— °- CDU)Cl) Cl)00DC•F C)t’)0(i-IC/) Cl) Cl) CD C) F- 0 D 0 -4,F CD CI)D C.— C D ciC/) D 0 D CDChacoCanyon200018501 65014501 2501050 85065020001850165014501250105085065001000020000300004000050000600007000000=0o.000—,0 -i-zrCD- =D0 C D CD D (0 0 D 0 :iJ9 Co C CD c)•WellsreportedInNaclmlento•WellsreportedInOjojJamoWellsreportedInQuatemary,8JIu’vulmChacoCanyon/•WellsreportedIn*WellsreportedIn+WellsreportedInKirtland/FruitlandFiciuredCllffsCliffHouseSanJuanRiver/CT22N124NT26NT28NOutcropofOjo,JamoSS•WellsreportedInNacimlento.WellsreportedInOjojejamoWellsreportedinQuaternaryAlluvuim•WellsreportedinKlrtland/FruhiandOutcropof*WellsreportedInFicturedCliffsActuredCliffsSS+WellsreportedinCliffHouse-1, CD C CD3’—0 0-Q).I_,—+aD03C,)-.Cl)0.-o C) CD CD()1D4CD a :3 CD U) 0 CD C,) :3 :3- CD U) a a U) U) C/) CD U) 0 :3 CO20001 8501650145012501050 8506502000185016501450125010508506500100002000030000400005000060000700000 0 D CD 0 D-n CD C -‘ CD 010’DC)0—•CDCDChacoCanyon/•WellsreportedkiNaclmlento•WellsreportedInOjoJamoSanJuanRiver/EEiE-j 0I122N••a124NaT26N128N-Il C CD a,a’— 0‘-JoCD_ D0 U)-.(fl-0 (00OD0C’) C) -‘ CD01CDOD CD ci D CD 0 -3 CD Cl) -4- 0 0 0 crOutcropofPicturedCliffsSSOutcropofOjoPJamoSS20001 8501 65014501 2501050 850650•WellsreportedInNaclmlento•WellsreportedInOjoJamo20001850165014501250105085065001000020000300004000050000600007000000a)=Figure 7.7 Conceptual cross section showing layering and model boundaryconditions.0G)00G)257-9CD C -‘ CDD1CD(D -4.J-‘QcTCDC,, CD a 0 Dr\):co3 CD CD CD Cf, -4- 3 CD U, C (I) CD CD Zr a -‘ C 0200018501650145012501050 850650200018501 650145012501050850650010000200003000040000500006000070000-1,Co C CD (0“—1(DD C) 0 Cl) -4- 0 0• (I) CDr\)I20001880200018501 650145012501050 8506501 650145012500105085010000200003000040000500006000070000650-9 C CD - 0 I CD ci 0 Cl) -‘ c5- C -4. 0 D cr cr CD 0 0 Dr\)0)Do-o-‘ 3 CD CD -‘ CD Cl) -4- 3 CD C,) -4, -‘ 0 3 0 D -4- 0 -I. (1) CD200018501 650145012501050 850650200018501 6501 45012501050850650010000200003000040000500006000070000a 0 0 0o LI) LI) LI)o aD CD-—0LI)(N0LI)00I-I)aD0If)(0000000aaCD0000It)00a00000F)0000(N00000o 0 0 0 0 0 0o LI) If) LI) LI) If) LI)0 CD (0 (N (0(N————Figure 7.11 The 14C distribution based on fina’ parameter estimates from joint dataset2610LI)aDCHAPTER 8. SUMMARYThe work in this thesis is devoted toward understanding parameter estimationfor groundwater flow models. Parameter estimation using only hydraulic head dataoften results in ill-conditioned parameter estimates, leading to problems with nonunique and unstable parameter estimates with large uncertainties. Prior informationon parameters or joint parameter estimation including tracer concentration data havebeen proposed to reduce the ill-conditioning of the parameter estimates. Priorinformation or joint estimation reduce the ill-conditioning of the parameter estimatesin some situations, but in other situation they do not reduce the ill-conditioning of theparameter estimates significantly. The focus of the thesis has been to understand inwhat situations prior information or joint estimation will significantly improve theparameter estimates, and develop guidelines for using prior information and jointparameter estimation in an efficient and responsible manner.The key contributions of this thesis are:1. Parameter space analysis using response surfaces and linearized confidenceregions is introduced. Response surfaces are used to understand conceptuallyhow prior information and joint parameter estimation improve parameterestimates. Linearized confidence regions are used to extend the concepts tomultiple parameter dimensions.2. Guidelines are developed for the efficient and responsible use of prior262information in groundwater flow models. Efficient use of prior informationinvolves identifying those parameters for which prior information will stabilizethe parameter estimates the most. Responsible use of prior informationinvolves those parameters for which errors in the prior information have thesmallest influence on the parameter estimates.3. Guidelines are developed to identify situations where joint parameter estimationwill significantly improve parameter estimates over single state parameterestimation. Criteria are also developed and evaluated for weighting individualdata sets in joint parameter estimation.4. The above guidelines are used to develop and calibrate a groundwater flowmodel for the San Juan Basin, New Mexico. Hydraulic head data, 14C data, andprior information on parameter values are used to calibrate the model. Themost efficient and responsible use of prior information is identified. Weightingcriteria are applied and evaluated to determine appropriate weights for the headand 14C data.263REFERENCESBaltz, E.H., Stratigraphy and regional tectonic implications of part of the UpperCretaceous and Tertiary rocks, east-central San Juan Basin, New Mexico, U.S.Geol. Surv. Prof. Pap., 552, 1967Baltz, E.H., S.R. Ash, and R.Y. Anderson, History of the nomenclature andstratigraphy of rocks adjacent to the Cretaceous-Tertiary boundary, westernSan Juan Basin, New Mexico, U.S. Geol. Surv. Prof. Pap., 524-D, 1966Campana, M.F., and E.S. Simpson, Groundwater residence times and discharge ratesusing discrete-state compartment model and 14C data, J. Hydrol., 72, 171-1 85,1984Carrera, J., State of the art of the inverse problem applied to the flow and solutetransport equations, in Groundwater Flow and Quality Modelling, edited by E.Custodio et al., pp. 549-583, D. Reidel Publishing Co., 1988Carrera, J., and S.P. Neuman, Estimation of aquifer parameters under transient andsteady state conditions, 1. Maximum likelihood method incorporating priorinformation, Water Resour. Res., 22(2), 199-210, 1986aCarrera, J., and S.P. Neuman, Estimation of aquifer parameters under transient andsteady state conditions, 2. Uniqueness, stability, and solution algorithms, WaterResour. Res., 22(2), 211-227, 1986bCarrera, J., and S.P. Neuman, Estimation of aquifer parameters under transient andsteady state conditions, 3. Aplication to synthetic and field data, Water Resour.Res., 22(2), 228-242, 1986cChavent, G., About the stability of the optimal control solution of inverse problems, inInverse and Improperly Posed Problems, edited by G. Anger, pp. 45-58,Akademie-Verlag, Berlin, 1979Cooley, R.L., A method of estimating parameters and assessing reliability of modelsof steady state groundwater flow, 1. Theory and numerical properties, WaterResour. Res., 13(2), 318-324, 1977Cooley, R.L., A method of estimating parameters and assessing reliability of modelsof steady state groundwater flow, 2. Application of statistical analysis, WaterResour. Res., 15(3), 603-617, 1979264Cooley, R.L., Incorporation of prior information on parameters into nonlinearregression groundwater flow models, 1. Theory, Water Resour. Res., 18(4),965-976, 1982Cooley, R.L., Exact Scheffe-type confidence intervals for output from groundwatermodels, 2. Combined use of hydrogeologic information and calibration data,Water Resour. Res., 29(1), 35-50, 1993Cooley, R.L. and R.L. Naff, Regression modeling of groundwater flow, USGS OpenFile Report 85-180, 1985Cooley, R.L., Konikow, L.F., and Naff, R.L., Nonlinear regression groundwater flowmodeling of adeep regional aquifersystem, Water Resour. Res., 22(13), 1759-1778, 1986Copty, N., Y. Rubin, and G. Mavko, Geophysical-hydrological identification of fieldpermeabilities through Bayesian updating, Water Resour. Res., 29(8), 2813-2825, 1993.Dagan, G., Stochastic modeling of groundwater flow by unconditional and conditionalprobabilities: The inverse problem, Water Resour. Res., 21(1), 65-72, 1985.Dagan, G., Time dependent macrodispersion for solute transport in anisotropicheterogeneous aquifers, Water Resour. Res., 24(9), 1491-1500, 1988Dagan, G., and Y. Rubin, Stochastic identification of recharge, transmissivity andstorativity in aquifer unsteady flow: A quasi-steady approach, Water Resour.Res., 24(10), 1698-1710, 1988.Emsellen, Y., and G. de Marsily, An automatic solution for the inverse problem, WaterResour. Res., 7(5), 1264-1 283, 1971Fasset, J.E., and J.S. Hinds, Geology and fuel resources of the Fruitland Formationand Kirtland Shale of the San Juan Basin, New Mexico and Colorado, U.S.Geol. Surv. Prof. Pap., 676, 1971Fontes, J. Ch., Environmental isotopes in groundwater hydrology, in Handbook ofEnvironmental Isotope Geochemistry, edited by P.Fritz and J. Ch. Fontes,pp.75-i 34, Elsevier, New York, 1980Fontes, J. Ch., and J.M. Gamier, Determination of the initial carbon-i 4 activity of thetotal dissolved carbon - A review of the existing models and a new approach,Water Resour. Res., 15, 399-41 3, 1979265Freeze, R.A., J. Massmann, L. Smith, T. Sperling, and B. James, Hydrogeologicaldecision analysis, 1. A framework, Ground Water, 28(5) 738-766, 1990Gailey, R.M., A.S. Crowe, and S.M. Gorelick, Coupled process parameter estimationand prediction uncertainty using hydraulic head and concentration data, Adv.Water Resour., 14(5), 301-314, 1991Gavalas,G.R., P.C. Shah, and J.H. Seinfeld, Reservoir history matching by Bayesianestimation, Soc. Pet. Eng., J., 16(6), 337-350, 1976Gelhar, L.W., and C.L. Axness, Three-dimensional stochastic analysis ofmacrodispersion in aquifers, Water Resour. Res., 19(1), 161 -180, 1983Graybill, F.A., Theory and Application of the Linear Model, Duxbury, North Scituate,Mass., 1976.Gupta, S.K., D. Lal, and P. Sharma, An approach to determining pathways andresidence time of groundwaters: Dual radiotracer dating, J. Geophys. Res., 86,5292-5300, 1981Hefez, E., V. Shamir, and J. Bear, Identifying the parameters of an aquifer cell model,Water Resour. Res., 11(6), 993-1004, 1975.Hill, M.C., A computer program (MODFLOWP) for estimating parameters of atransient, three dimensional ground-water flow model using nonlinearregression, U.S. Geol. Sur. Open File Rep. 91 -484, Denver, 1992Hill, P.D., A review of experimental design procedures for regression and modeldiscrimination, Technometrics, 20, 15-21,1978Hoeksema, R.J., and P.K. Kitanidis, An application of the geostatistical approach tothe inverse problem in two-dimensional groundwater modeling, Water Resour.Res., 20(7), 1003-1 020, 1984Jacquard, P., and C. Jam, Permeability distribution from field pressure data, Soc. Pet.Eng. J., 5(4), 281-294, 1965.Kitanidis, P.K., and E.G. Vomvoris, A geostatistical approach to the inverse problemin groundwater modeling and one-dimensional simulations, Water Resour. Res.,19(3), 677-690, 1983.Kleinecke, D., Use of linear programming for estimating geohydrologic parameters ingroundwater basins, Water Resour. Res., 7(2), 367-375, 1971.266Krabbenhoft, D.P., M.P. Anderson, C.J. Bowser, and J.W. Valley, Estimatinggroundwater exchange with lakes, 1. The stable isotope mass balance method,Water Resour. Res. 26(10), 2445-2453, 1990a.Krabbenhoft, D.P., M.P. Anderson, and C.J. Bowser, Estimating groundwaterexchange with lakes, 2. Calibration of a three-dimensional, solute transportmodel to a stable isotope plume, Water Resour. Res. 26(10), 2455-2462,1990b.Kuczera, G. Improved parameter inference in catchment models, 1. Evaluatingparameter uncertainty, Water Resour. Res., 19(5), 1151 -1162, 1983a.Kuczera, G. Improved parameter inference in catchment models, 2. Combiningdifferent types of hydrologic data and testing their compatability, Water Resour.Res., 19(5), 1163-1172, 1983b.Marquardt, D.W., An algorithm for least-squares estimation of non-linear parameters,J. Soc. md. Appi. Math., 11(2), 431-441, 1963Mendoza, C.A., R. Therrien, and E.A. Sudicky, ORTHOFEM users guide, Universityof Waterloo, Waterloo, Ontario, l4pp.,1992Mishra, S., and S.C. Parker, Parameter estimation for coupled unsaturated flow andtransport, Water Resour. Res., 25(3), 385-396, 1989.Mook, W.G., Carbon-14 in hydrogeological studies, in Handbook of EnvironmentalIsotope Geochemistry, edited by P.Fritz and J. Ch. Fontes, pp.49-74, Elsevier,New York, 1980Neuman, S.P. and S. Yakowitz, A statistical approach to the inverse problem inaquifer hydrology, 1. Theory, Water Resour. Res., 15(4), 845-860, 1979.Neuman, S.P., A statistical approach to the inverse problem of aquifer hydrology, 3.Improved solution method and added perspective, Water Resour. Res., 16(2),331-346, 1980.Nour-Omid, B., W.S. Dunbar, and A.D. Woodbury, Lanczos and Arnoldi methods forthe solution of convection-diffusion equations, Comp. Meth. App. Mech. Eng.,88, 75-95, 1991Pearson, F.J., Jr., and D.E. White, Carbon-14 ages and flow rates of water in theCarrizo Sand, Atascosa County, Texas, Water Resour. Res., 3, 251-261, 1967267Phillips, F.M., M.K. Tansey, L.A. Peeters, S. Cheng, and A. Long, An isotopicinvestigation of groundwater in the central San Juan basin, New Mexico:Carbon-14 dating as a basis for numerical flow modeling, Water Resour. Res.25(10), 2259-2273, 1989Phillips, F.M. and M.K. Tansey, An integrated isotopic/physical approach to anumerical model of groundwater flow in the San Juan Basin, Rep. 197, NewMexico Water Resour. Res. Inst., Las Cruces, 1984Press, W.H., B.P.Flannery, S.A. Teukolsky, and W.T. Vettering, Numerical Recipes,Cambridge University Press, Cambridge, 1986Ratowsky, D.A., Nonlinear Regression Modeling, Marcel Dekker, New York. 1984Robertson, W.D. and J. Cherry, Tritium as an indicator of recharge and dispersion ina groundwater system in central Ontario, Water Resour. Res. 25(6),1097-1109, 1989Shah, P.C., G.R. Gavalas, and J.H. Seinfeld, Error analysis in history matching: Theoptimum level of parameterization, Soc. Pet. Eng. J., 18(3), 21 9-228, 1978Sorooshian, S., and F. Arfi, Response surface sensitivity analysis methods for postcalibration studies, Water Resour. Res., 18(5), 1531 -1 538, 1982.Sorooshian, S., and V.K. Gupta, Automatic calibration of conceptual rainfall-runoffmodel parameters: The question of parameter observability and uniqueness,Water Resour. Res., 19(1), 260-268, 1983.Sorooshian, S., and V.K. Gupta, The analysis of structural identifiability: Theory andapplication to conceptual rainfall-runoff models, Water Resour. Res., 21(4),487-495, 1985.Stallman, R.W. Numerical analysis of regional water levels to define aquifer hydrology,Trans. AGU, 37(4), 451-460, 1956Stone, W.J., F.P. Lyford, P.F. Frenzel, N.H. Mizell, and E.T. Padgett, Hydrogeologyand water resources of the San Juan Basin, New Mexico, Miner. Resour. Rep.N.M. Bur. Mines Miner. Resour.,6, 7Opp,1983Strecker, E.W., and W. Chu, Parameter identification of a groundwater contaminanttransport model, Ground Water, 24, 56-62, 1986.Sun, N., and W.W-G. Yeh, Identification of parameter structure in groundwater inverseproblems, Water Resour. Res., 21(6), 869-883, 1985268Sun, N., and W.W-G. Yeh, Coupled inverse problems in groundwater modeling, 1.Sensitivity analysis and parameter identification, Water Resour. Res., 26(10),2507-2525, 1990aSun, N., and W.W-G. Yeh, Coupled inverse problems in groundwater modeling, 2.Identifiability and experimental design, Water Resour. Res., 26(10), 2526-2551,1990bTarantola, A. and B. Valette, Generalized nonlinear inverse problems using the leastsquares criterion, Rev. Geophys., 20(2), 21 9-232, 1982Toorman, A.F., P.J. Wierenga, and R.G. Hills, Parameter estimation of hydraulicproperties from one step outflow data, Water Resour. Res., 28(11), 3021-3028,1992Umari, A.R., R. Willis, and P.L.-F. Liu, Identification of aquifer dispersivities in twodimensional transient groundwater contaminant transport: An optimizationapproach, Water Resour. Res., 15(4), 81 5-831, 1979.Vecchia, A.V., and R.L. Cooley, Simultaneous confidence and prediction intervals fornonlinear regression models with application to a groundwater flow model,Water Resour. Res., 23(7), 1237-1 250, 1987.Vogel, J.C., Investigation of groundwater flow with radiocarbon, in Isotopes inHydrology, pp. 235-237, International Atomic Energy Agency, Vienna, 1967Wagner, B.J. and S.M. Gorelick, A statistical methodology for estimating transportparameters: Theory and applications to one dimensional advective-dispersivesystems, Water Resour. Res., 22(8), 1303-131 5, 1986.Wagner, B.J. and S.M. Gorelick, Optimal groundwater quality management underparameter uncertainty, Water Resour. Res., 23(7), 1162-1174, 1987.Wang K., P.Y. Shen, and A.E. Beck, A solution to the inverse problem of coupledhydrological and thermal regimes, in Hydrological regimes and their subsurfacethermal effects, edited by A.E. Beck, L. Stegena and G. Garven, AGUMonograph series, Washington, 1989Weiss, R. and L. Smith, Parameter estimation using hydraulic head and environmentaltracer data, in Proceedings, 1993 Ground Water Modeling Conference, Golden,Colorado, 60-70, 1993269Woodbury, A.D., and L. Smith, Simultaneous inversion of hydrogeologic and thermaldata, 2. Incorporation of thermal data, Water Resour. Res., 24(3), 356-372,1988Woodbury, A.D., L. Smith, and W.S. Dunbar, Simultaneous inversion of hydrogeologicand thermal data, 1. Theory and application using hydraulic head data, WaterResour. Res., 23 (8), 1586-1 606, 1987Woodbury, A.D., W.S. Dunbar, B. Nour-Omid, Application of the Arnoldi algorithm tothe solution of the àdvection-dispersion equation, Water Resour. Res., 26(10),2579-2590, 1990 -Xiang, Y., J.F. Sykes, and N.R. Thomson, A composite Li parameter estimator formodel fitting in groundwaterflow and solute transport simulation, Water Resour.Res., 29(6), 1661-1673, 1993Yakowitz, S. and L. Duckstein, Instability in aquifer identification: Theory and casestudies, Water Resour. Res., 16(6), 1045-1064, 1980Yeh, W. W-G., Optimal identification of parameters in an inhomogeneous medium withquadratic programming, Soc. Pet. Eng. J., 15(5), 371 -375, 1975.Yeh, W. W-G., and Y.S. Yoon, Parameter identification with optimum dimension inparameterization, Water Resour. Res., 17(3), 664-672, 1981Yeh, W. W-G., Y.S. Yoon, and K.S. Lee, Aquifer parameter identification with krigingand optimum parameterization, Water Resour. Res., 19(1), 225-233, 1983.Yeh, W. W-G., Review of parameter identification procedures in groundwaterhydrology: The inverse problem, Water Resour. Res., 22 (2), 95-1 08, 1986270
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Efficient and responsible use of prior information...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Efficient and responsible use of prior information and joint parameter estimation for estimating parameters… Weiss, Richard A. 1994
pdf
Page Metadata
Item Metadata
Title | Efficient and responsible use of prior information and joint parameter estimation for estimating parameters of groundwater flow models |
Creator |
Weiss, Richard A. |
Date Issued | 1994 |
Description | Parameter estimation problems for groundwater flow systems are often non-identifiable or unstable using only hydraulic head data. Prior information on parameters or joint parameter estimation including tracer concentrations may be used to stabilize the parameter set. Response surfaces and multiparameter confidence regions are used to identify the most efficient and responsible use of prior information and joint data sets for the purpose of parameter estimation. Prior information on some of the parameters are used to stabilize the parameter set for parameter estimation. Efficient use of prior information involves identifying those parameters for which prior information will stabilize the model parameter set the most. Responsible use of prior information involves identifying how errors in the prior information will influence the parameter estimates. The most responsible parameters for prior information are those parameters for which errors in the prior information have the least influence on the final parameter estimates. Guidelines are developed for efficient and responsible use of prior information in parameter estimation based on an analysis of the parameter space using response surfaces and eigenspace decomposition. Joint parameter estimation is used when more than one data set is available to estimate the model parameters. Response surfaces and confidence regions are used to show how multiple data sets reduce parameter uncertainty. The value of future data in reducing the uncertainty of parameter estimates is explored. For weighting data sets in joint parameter estimation, three criteria based on parameter space analysis are proposed. These three criteria are evaluated and compared to the more traditional weighting method based on an analysis of data residuals. A groundwater model for the San Juan Basin, New Mexico, is constructed and calibrated using the methods developed in this thesis. Hydraulic head data, ¹⁴C concentration data, and prior information on the model parameters is used to calibrate the model in an efficient and responsible manner. The model is calibrated in two stages. In the first stage, prior information on the hydraulic conductivity parameters for the lower model layers was found to be both efficient and responsible in stabilizing the parameter set. In the second stage, the parameter estimates and uncertainties based on the four weighting criteria were compared. |
Extent | 5290959 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-04-15 |
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. |
IsShownAt | 10.14288/1.0052915 |
URI | http://hdl.handle.net/2429/7178 |
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 |
GraduationDate | 1994-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1994-954086.pdf [ 5.05MB ]
- Metadata
- JSON: 831-1.0052915.json
- JSON-LD: 831-1.0052915-ld.json
- RDF/XML (Pretty): 831-1.0052915-rdf.xml
- RDF/JSON: 831-1.0052915-rdf.json
- Turtle: 831-1.0052915-turtle.txt
- N-Triples: 831-1.0052915-rdf-ntriples.txt
- Original Record: 831-1.0052915-source.json
- Full Text
- 831-1.0052915-fulltext.txt
- Citation
- 831-1.0052915.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:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052915/manifest