USGS Publications Warehouse

Abstract

Construction of a ground-water model for a field area is not a straightforward process. Data are virtually never complete or detailed enough to allow substitution into the model equations and direct computation of the results of interest. Formal model calibration through optimization, statistical, and geostatistical methods is being applied to an increasing extent to deal with this problem and provide for quantitative evaluation and uncertainty analysis of the model. However, these approaches are hampered by two pervasive problems: 1) nonlinearity of the solution of the model equations with respect to some of the model (or hydrogeologic) input variables (termed in this report system characteristics) and 2) detailed and generally unknown spatial variability (heterogeneity) of some of the system characteristics such as log hydraulic conductivity, specific storage, recharge and discharge, and boundary conditions. A theory is developed in this report to address these problems. The theory allows construction and analysis of a ground-water model of flow (and, by extension, transport) in heterogeneous media using a small number of lumped or smoothed system characteristics (termed parameters). The theory fully addresses both nonlinearity and heterogeneity in such a way that the parameters are not assumed to be effective values. The ground-water flow system is assumed to be adequately characterized by a set of spatially and temporally distributed discrete values, ?, of the system characteristics. This set contains both small-scale variability that cannot be described in a model and large-scale variability that can. The spatial and temporal variability in ? are accounted for by imagining ? to be generated by a stochastic process wherein ? is normally distributed, although normality is not essential. Because ? has too large a dimension to be estimated using the data normally available, for modeling purposes ? is replaced by a smoothed or lumped approximation y?. (where y is a spatial and temporal interpolation matrix). Set y?. has the same form as the expected value of ?, y 'line' ? , where 'line' ? is the set of drift parameters of the stochastic process; ?. is a best-fit vector to ?. A model function f(?), such as a computed hydraulic head or flux, is assumed to accurately represent an actual field quantity, but the same function written using y?., f(y?.), contains error from lumping or smoothing of ? using y?.. Thus, the replacement of ? by y?. yields nonzero mean model errors of the form E(f(?)-f(y?.)) throughout the model and covariances between model errors at points throughout the model. These nonzero means and covariances are evaluated through third and fifth-order accuracy, respectively, using Taylor series expansions. They can have a significant effect on construction and interpretation of a model that is calibrated by estimating ?.. Vector ?.. is estimated as 'hat' ? using weighted nonlinear least squares techniques to fit a set of model functions f(y'hat' ?) to a. corresponding set of observations of f(?), Y. These observations are assumed to be corrupted by zero-mean, normally distributed observation errors, although, as for ?, normality is not essential. An analytical approximation of the nonlinear least squares solution is obtained using Taylor series expansions and perturbation techniques that assume model and observation errors to be small. This solution is used to evaluate biases and other results to second-order accuracy in the errors. The correct weight matrix to use in the analysis is shown to be the inverse of the second-moment matrix E(Y-f(y?.))(Y-f(y?.))', but the weight matrix is assumed to be arbitrary in most developments. The best diagonal approximation is the inverse of the matrix of diagonal elements of E(Y-f(y?.))(Y-f(y?.))', and a method of estimating this diagonal matrix when it is unknown is developed using a special objective function to compute 'hat' ?. When considered to be an estimate of f(y?.), the estimate f(y 'hat' ?) is biased because of nonlinearity in f(?) and f(y?) (?= ?. or 'hat' ?), but when considered to be an estimate of f(?), f(y 'hat' ?) is biased only because of components of nonlinearity in f(?) and f(y?) known as the intrinsic nonlinearity. (Intrinsic nonlinearity in either f(?) or f(y?) is a component of the total nonlinearity in either function that cannot be eliminated by a unique transformation of either (? or ?, respectively.) Because both types of intrinsic nonlinearity can be small, f(y 'hat' ?) can be nearly unbiased as an estimate of f(?) . Analogous results hold for a prediction g(y 'hat' ?) (where g is some function of parameters of interest to the investigator), except that in this case the intrinsic nonlinearity is for the combination of either g(?) and f(?) or g(y?) and f(y?) , termed the combined intrinsic nonlinearity. The biases are evaluated to second-order accuracy using Taylor series expansions and the analytical least squares solution, but an investigator would probably be more interested in estimates of f(?) and g(?) than in estimates of fictitious variables f(y?.) and g(y?.). Predictive accuracy of a model is thus strongly tied to the degree of intrinsic nonlinearity of the models f(?) and f(y?.) together with the combined intrinsic nonlinearity of the models and the predictions to be made with them. Uncertainties in the estimates of ?., f(y?.), f(?), g(y?.), and g(?) (or some future measurement of g(?) are addressed through nonlinear confidence regions, confidence intervals, and prediction intervals. If ? and the observation errors are normally distributed, statistical distributions of functions of the weighted sums of squared errors in the estimates necessary to define the regions and intervals approximate F distributions that are modified with correction factors. These functions are very similar to the standard ones developed for linear models except that the parameter set ?. is stochastic rather than fixed. The correction factors correct the distributions to account for intrinsic nonlinearity of f(y?) and deviation of the weight matrix from the correct one. The correction factors are derived using the Taylor series and perturbation method used for the least squares; the generality of the factors and concepts leading to them are verified by an independent method that does not rely on Taylor series and perturbations. Because of the effects of spatial correlation, confidence regions and confidence intervals would generally be too small without using components of the correction factors needed to correct for using an incorrect weight matrix unless the correct one is used; prediction intervals may often be nearly correct. Approximate bounds ,for the correction factors are developed for use when the information on nonlinearity and heterogeneity necessary to calculate them is not available. Measures of total model nonlinearity of f(y?), intrinsic nonlinearity of f(y?), and combined intrinsic nonlinearity of (y?) and g(y?) help an investigator decide when the components of the correction factors accounting for the types of intrinsic nonlinearity are not important. Two examples are analyzed to test the validity and robustness of the theory when the model error is large. Example 1 is for one-dimensional, steady-state flow in an aquifer having log transmissivity (lnT) that varies stochastically at small scale and recharge (W) that is constant. Example 2 is for two-dimensional, steady-state flow in a zoned aquifer where lnT and W vary spatially at both large and small scales, the small-scale variations being stochastic. Hydraulic head data Y were generated as f(y?.) plus the sum of model errors f(?)-f(y?.) and small, zero-mean, independent, normal observation errors, Y-f(?). The most important results are as follows: 1) The total nonlinearity in (y?) is large for both examples, but the intrinsic nonlinearity in f(y?) and the combined intrinsic nonlinearity of f(y?) and g(y?) for g(y?) equal to lnT, W, and a predicted hydraulic head, are all small for both examples. As the theory predicts, the corresponding biases also were found to be small. 2) The sum of the model and observation errors, Y-f(y?.), has a nonzero mean and is not normally distributed for either example. 3) Spatial correlations between dements of Y-f(y?.) are often large, and large values are usually positive for both examples. 4) Residual set Y-f(y 'hat' ?) does not differ significantly from the zero-mean, normally distributed set predicted by the theory for either example. 5) For example 2, use of the correct, full weight matrix produced accurate confidence and prediction intervals, as determined using a Monte Carlo procedure. (Because of severe ill conditioning, example 1 produced open-ended intervals that could not be formally analyzed.) 6) Use of the best diagonal matrix for example 2 produced accurate confidence intervals only when the correction factors were used. Otherwise, the intervals are too small as determined by the Monte Carlo procedure. The prediction intervals did not have to be corrected. 7) Use of the estimated weight matrix for example 2 produced confidence and prediction intervals that had to be corrected, but the minimum containment probability of 0.92 for the intervals after correction remains slightly too small compared to the nominal probability of 0.95 because the correction factors are more approximate for this case than when the best diagonal weight matrix is used.