Unrealistic parameter estimates in inverse modelling - Hydrologie.org

0 downloads 0 Views 807KB Size Report
Abstract Estimation of unrealistic parameter values by inverse modelling ... the true conditions and character of the errors are completely known. ..... McKenna, S. A. & Poeter, E. P. (1995) Field example of data fusion for site characterization.
Calibration and Reliability in Groundwater Modelling (Proceedings of the ModelCARE 96 Conference held at Golden, Colorado, September 1996). IAHS Publ. no. 237, 1996.

277

Unrealistic parameter estimates in inverse modelling: a problem or a benefit for model calibration? EILEEN P. POETER Department of Geology and Geological Engineering, Colorado School of Mines, 1500 Illinois Street, Golden, Colorado 80401, USA

MARY C. HILL US Geological Survey, Water Resource Division, Box 25046, MS 413, Lakewood, Colorado 80225, USA

Abstract Estimation of unrealistic parameter values by inverse modelling is useful for constructed model discrimination. This utility is demonstrated using the three-dimensional, groundwater flow inverse model MODFLOWP to estimate parameters in a simple synthetic model where the true conditions and character of the errors are completely known. When a poorly constructed model is used, unreasonable parameter values are obtained even when using error free observations and true initial parameter values. This apparent problem is actually a benefit because it differentiates accurately and inaccurately constructed models. The problems seem obvious for a synthetic problem in which the truth is known, but are obscure when working with field data. Situations in which unrealistic parameter estimates indicate constructed model problems are illustrated in applications of inverse modelling to three field sites and to complex synthetic test cases in which it is shown that prediction accuracy also suffers when constructed models are inaccurate.

INTRODUCTION Developing conceptual models, choosing a computer code and developing model designs are the preliminary steps in a groundwater modelling project after first establishing the purpose of the project (Anderson & Woessner, 1992). These steps define the numerical models representing a given field situation, including: (a) the most important sources and sinks of water in the field system and how they are to be simulated; (b) the available data on the geohydrological system; (c) the system geometry (generally the number and type of model layers and the areal extent of these layers); (d) the spatial and temporal structure of the hydraulic properties (generally using zones of constant value or deterministic or stochastic interpolation methods); and (e) boundary condition location and type. We call the resulting models constructed models. To produce acceptable models, the best set of parameter values (and associated confidence intervals) of each constructed model need to be determined. In addition, the best calibrated models need to be identified. Thus, there are three levels of model calibration: (a) model construction; (b) parameter estimation; and (c) model discrimination. While for some constructed models the parameter estimation is unique, the overall model calibration is never unique for the complicated groundwater systems commonly

278

Eileen P. Poeter & Mary C. Hill

considered. In practice, multiple constructed models must be developed at the beginning of a project and, depending on the character of the incoming data and results of ongoing analyses, each model is either retained for further consideration or eliminated during the modelling process. Parameter estimation can be approached using an inverse model (Poeter & Hill, 1996); for example, nonlinear regression can be used to find the set of parameter values that provides the best fit of model results to field observations, where "best fit" is defined as minimizing the value of the sum-of-squared weighted residuals. However, parameter estimation is often accomplished by a trial-and-error approach, during which the modeller iteratively selects parameter values to improve the fit of model results to field observations using intuition about model response to changes in parameters, and knowledge of reasonable parameter ranges. The time consuming nature of intuitive parameter value adjustment limits the range of alternative constructed models that are considered and, given the lack of rigorous analysis of parameter correlations, variance/ covariance, and residuals, there is no assurance that the estimated parameter values for any model are "the best". Consequently, conclusive model discrimination is nearly impossible. This shortcoming is often so extreme that only one constructed model is considered. With inverse models used to determine parameter values that optimize the fit of the model results to the field observations for a given model configuration, the modeller is freed from tedious trial-and-error calibration involving changes in parameter values so more time can be spent addressing insightful questions about the hydrological system. If a constructed model is not an adequate representation of the groundwater flow system, then the estimated parameter values are not likely to reflect field conditions, predictions made using the model are likely to be in error, and decisions based on those predictions may not produce the best, or even a reasonably accurate, result. Indicators of inaccurately constructed models also include non-random weighted residuals and poor fit to the data, but the indicator focused on in this work is unreasonable parameter values (e.g. a lower hydraulic conductivity for a sand than a clay, or a boundary flux with the wrong sign) yielding the best fit between the observed and simulated values (such as hydraulic heads, concentrations and flows). Thus, during calibration it is important to determine if the best model fit is achieved for unreasonable parameter values and whether these parameter values are well estimated with the available data. Best fitting unreasonable parameter values are more likely to be discovered using the inverse modelling approach than by trial and error, because the inverse model will determine the parameters that provide the best fit evei though they may be of unreasonable absolute or relative magnitude. A modeller using the trial-and-error approach generally will not try unreasonable values (e.g. a much higher hydraulic conductivity for a silt deposit than an adjacent sandy-gravel deposit). Instead, usually unknowingly, the modeller accepts greater discrepancies in the match between field observations and model results (i.e. sacrifices the best fit between the model and the data) to maintain reasonable model parameters. The output of unreasonable values is often disillusioning to new users of inverse models, but it is actually a great benefit. At first glance it may seem that the trial-anderror approach provides the more reasonable result and so would be the preferred method. However, when the modeller sacrifices honouring data to maintain "reasonableness", important information about the system is ignored and evidence of error in the constructed model is not allowed to surface. It is our experience that well estimated,

Unrealistic parameter estimates in inverse modelling

279

unreasonable parameter values produced by the inverse model indicate that the constructed model is incorrect. With this information and the statistics calculated by the inverse model, the modeller can call upon experience with the geohydrological systems and knowledge of reasonable parameter values to develop and test numerous constructed models, increasing the likelihood that accurate models will be used for prediction.

IMPACT OF THE CONSTRUCTED MODEL ON PARAMETER ESTIMATION Synthetic problems are useful for illustrating the impact of a poorly constructed model on the estimated parameter values because the "true" subsurface conditions are known and can be compared to the outcome of the modelling process.

Synthetic model configuration A simple basin configuration, roughly 30 X 40 km is presented in Fig. 1(a). Flow is toward the streams and northward with groundwater discharge to streams and a body of surface water occupying three cells at the north end of the basin. The perimeter of the basin is defined by the groundwater divide in the east, west, and south. The saturated (a) Hydraulic Head Distribution

•5000

0

5000

10000

Fig. 1 Simple deterministic synthetic basin configuration: (a) hydraulic head distribution; (b) hydraulic conductivity (0.1, to 1, 5 and 10 m day"1 respectively from black to lightest grey); and (c) recharge increases linearly with decreasing grey tone from 1.5 x 10"5 to 3 X 10'4 m day"1, black is zero recharge or outside the basin boundary.

280

Eileen P. Poeter & Mary C. Hill

Observation Locations

b

Base Case (Cases 1 through 4)

c

Case 5

Fig. 2 Observation locations and alternative zonations for hydraulic conductivity.

thickness of the unconfined aquifer varies from 30 to 150 m within the basin. This synthetic basin is hydraulically two-dimensional: that is, there is no variation of hydraulic character with depth. Consequently, the hydraulic conductivity distribution (analogous to hydrofacies distribution) shown on the map prevails to the depth at which bedrock occurs. Four hydrofacies are delineated with increasing mean grain size, and hydraulic conductivity (K) increases from 0.1, to 1, 5 and 10 m day"1 respectively from black to lightest grey in Fig. 1(b). Recharge is in the form of infiltration, varying from 1.5 X 10"5 m day"1 in the south to 3 X 10"4 m day"1 in the northern area (Fig. 1(c)). Recharge is zero in the regional discharge area at the north end of the basin where vegetative consumption prevents infiltration to the water table. Groundwater flow is simulated in the synthetic basin to provide true values of hydraulic heads (Fig. 1(a)). True hydraulic head and hydraulic conductivity at 16 observation locations (Fig. 2(a)) and one true groundwater discharge measurement of 0.2 m3 s"1 to the northeastern tributary are used in the inverse modelling illustration to estimate hydraulic conductivity throughout the area and recharge over the basin. Parameter estimation using alternative constructed models We use MODFLOWP (Hill, 1992) to estimate parameters and find it to be robust. In

Unrealistic parameter estimates in inverse modelling

281

calibrating the simple basin model to the observed heads and flows, we find that errors in field observations, in the initial estimates of parameters or in minor variations in zonation of geohydrological units do not cause problems in obtaining a reasonable parameter estimation solution. Problems arise when a poorly constructed model (e.g. zonation significantly different from the true zonation) is hypothesized. In such a case, unreasonable parameter values are obtained even when error free observations and true starting values are used for the parameters. Although this may be perceived as a problem, it is actually beneficial because it renders inverse modelling an excellent tool for differentiating accurate and inaccurate constructed models. For this simple basin model, the hydraulic conductivity value of the four hydrofacies (K1-K4) and the magnitude of the recharge R (given its relative spatial distribution) are estimated. Hydraulic conductivity distributions and estimated parameter values for the following cases are presented in Figs 2 and 3 respectively. Case 1 - the base case: with error free observations, a perfect constructed model, and correct "true" values of parameters as a starting point, the parameter estimation process quickly converges on the true parameter values (Figs 2(b) and 3). Case 2 - incorrect initial estimates: with error free observations, a perfect constructed model and incorrect initial values of the parameters including various combinations of parameters defined as two of orders of magnitude too small or large, the parameter estimation process quickly converges on the true parameter values (Figs 2(b) and 3). Case 3 - error in observations: with error laden observations having normally distributed noise with a standard deviation of 2 m on the head observations and the flow underestimated by 4%, a perfect constructed model and true parameter values for initial estimates, the code converges readily to values near the true values. The estimated parameter values are slightly different from the true values, reflecting parameter values that better fit the erroneous observations (Figs 2(b) and 3).

Parameter Estimates

7 w

6

n E 5

z

a> 4 tn to

"

3 2 1 n 0.001

0.01 0.1 Parameter Value (m/day)

1 10 [recharge x 100]

100

Fig. 3 Estimated parameter values for alternative models presented in Fig. 2.

282

Eileen P. Poeter & Mary C. Hill

Case 4 - larger error in observations: identical to Case 3, but with a standard deviation of 4 m on head observations and a flow measurement underestimated by 7.5 % (Figs 2(b) and 3). Case 5 - minor zonation error no. 1: with error free observations and accurate initial values for the parameters, a more discontinuous definition of the coarsest grained (high hydraulic conductivity) zone than exists in the base case (Fig. 2(c)) results in rapid convergence on parameter values close to the correct values (Fig. 3). Case 6 - minor zonation error no. 2: when the percentage of finest grained facies is overestimated and the coarsest grained facies is connected to the hydraulic boundary on the north end of the domain (Fig. 2(d)), errors in parameter estimations begin to arise. In order to match the observed heads better, hydraulic conductivity of the fine grained facies (Kl) is overestimated, compensating for the overabundance of the facies in the model. Hydraulic conductivity is underestimated for the coarse grained facies (K4), compensating for the hydraulic connection that does not actually exist (Fig. 3). The cases presented to this point pose some formidable problems to the parameter estimation algorithm without difficulty in obtaining reasonable estimates. The remaining two cases show how larger errors in the facies zonation pattern can create large errors in the estimates of parameter values. Case 7 - major zonation error: when the percentage of fine grained facies (Kl) is substantially over estimated and the percentage of coarse grained facies (K4) is underestimated and is much more discontinuous than in the base case (Fig. 2(e)), an unsatisfactory set of parameter estimates results. All of the parameters, including recharge rate, are estimated to be nearly an order of magnitude or more lower than the base case values. The estimated K of the coarse grained facies is lower than the medium grained facies (Fig. 3). Case 8 — major zonation error: sometimes the predominant geology in a small borehole is not thought to represent the geology of a 4 x 106 m2 flow model grid block, so this variation includes three locations where the geological observations are not honoured for the flow model grid block (Fig. 2(f)) and yields a more reasonable recharge rate, but overestimates the K of the three finest grained facies. The relative order of hydraulic conductivity between facies is again incorrect with the finest grained facies exhibiting a higher hydraulic conductivity than the coarsest grained facies and a reversal in the expected Ks for the two coarsest grained facies (Figs 3).

EXAMPLES OF MODEL CONSTRUCTION PROBLEMS IN COMPLEX SYSTEMS The methods described above for a simple test case are powerful in the analysis of more complex systems. This has been demonstrated using complex synthetic test cases by Poeter & McKenna (1995) and in the analysis of field sites using groundwater flow

Unrealistic parameter estimates in inverse modelling

283

models of Otis Air Force Base (Anderman et al., 1996) and a Colorado School of Mines Test Site (McKenna & Poeter, 1995), and using a convective-dispersive model of the Grindsted Old Landfill in Denmark (Christiansen et al, 1995; Barbelo et al, 1996). Highlights from these analyses are summarized below. Complex synthetic test cases The three-dimensional test cases investigated by Poeter & McKenna (1995) using MODFLOWP demonstrate that the conclusions drawn above for the simple test case are also valid for a more complicated system characterized by three-dimensional flow and substantial heterogeneity. This further supports the value of application of inverse modelling to field problems. Otis Air Force Base The two-dimensional areal model of LeBlanc (1984), which was calibrated by trial and error, was recalibrated with nonlinear regression methods to test a new method of using concentration data in inverse modelling (Anderman et al, 1996). A version of MODFLOWP modified to include observations of the path and time of advective travel was used. The observed path and time of advective travel was inferred from concentration measurements of a sewage plume that was introduced into the groundwater system in the early 1940s. Other data included in the regression were hydraulic heads and net flow into the groundwater system from Ashumet pond. Estimated parameters included a homogeneous aquifer hydraulic conductivity, the hydraulic conductivity of the bed of Ashumet Pond, flow rates along the northern model boundary and at the sewage disposal site, and a uniform areal recharge rate. Using the model construction described by LeBlanc (1984), it was found that the best fit to the advective travel, head and flow data was achieved with reasonable parameter values except for the areal recharge rate, which was about half the expected rate. In addition, the recharge rate was estimated precisely enough that a linear 95% confidence limit interval constructed about the estimate did not even come close to including reasonable values. This indicates that the data used are sufficient to distinguish between an accurate and inaccurate model, and that the constructed model was significantly inaccurate in some way. A variety of potential model construction problems that might cause the unrealistic recharge rate were tested by using regression to find the best fit parameter values in each case. One of the considered alternatives proved to be a plausible explanation of the problem; it involved the constant head boundary condition imposed along the southern boundary and southern parts of the east and west boundaries of the model. These boundaries represent surface water bodies. The elevations of the surface water bodies were derived from 10 foot contour topographic maps, and these elevations were used as the defined head at these boundaries, as is common practice in the development of groundwater models. It was found that if these "measured" heads were consistently higher (by just 1 foot) than the surface water body levels when the hydraulic heads used in the regression were measured, the estimated recharge rate would be about half of a reasonable rate. This indicates that, in systems such as this in which

284

Eileen P. Poeter & Mary C. Hill

surface water bodies define boundary conditions through which most or all of the water flows, and groundwater levels are only on the order of 5 feet above the surface water bodies, it is important to measure the elevation of the surface water bodies accurately at the same time as groundwater heads are measured. Colorado School of Mines test site McKenna & Poeter (1995) combined geological and geophysical information into a description of hydrological heterogeneity. The relationship between hydrofacies (defined using geological description of cuttings and cores, geophysical logs and permeability measurements of core and from packer tests) and seismic velocity was defined subjectively and discriminant analysis was used to determine the probability that a given location belonged to the hydrofacies. The hydrofacies were described with indicator values. Multiple indicator, geostatistical realizations of hydrofacies at locations between boreholes were generated first using only hard data (data with negligible uncertainty, i. e. those with greater than 95% probability of belonging to the hydrofacies), and again, supplementing these data with soft data (data with non-negligible uncertainty, i.e. those with less than 95 % probability of belonging to a hydrofacies and seismic tomography measurements). The plausibility of each stochastic hydrofacies zonation was determined via parameter estimation using MODFLOWP. Use of soft data, coupled with elimination of realizations when parameter estimation revealed a poor fit and/or unreasonable parameter values, resulted in narrower confidence limits on the estimated values. This sensitivity to fine scale geologically based zonation patterns is fortunate because it allows use of hydraulic head data and prior information on reasonable parameter values to delineate the small scale heterogeneity that is critical to the migration of contaminants and reduces the uncertainty associated with predicted flow through the site. Grindsted Old Landfill This study is the first to use nonlinear regression in the calibration of a threedimensional advective-transport model (Christiansen et al, 1995). The system is characterized by a steady state flow field and layered, apparently homogenous hydrographic units. Data included in the regression were hydraulic heads and concentrations. Estimated parameters included horizontal and vertical hydraulic conductivities for three relatively permeable layers, the vertical hydraulic conductivity of a confining unit, and the longitudinal dispersivity. Regression results using only the hydraulic head data produced unreasonable estimates of several parameters, but the confidence intervals on these parameter estimates were extremely large and included reasonable values. This indicated that the head data were insufficient to determine whether model construction error was a problem. When concentration data were also included in the regression, the residuals were unbiased, normally distributed, and were randomly distributed in space, and all estimated parameter values were reasonable, indicating that model construction was reasonably accurate.

Unrealistic parameter estimates in inverse modelling

285

REFERENCES Anderman, E. R., Hill, M. C. & Poeter.E. P. (1996) Two-dimensional advective transport in groundwater flow parameter estimation. Groundwater (in press). Anderson,M. P. &Woessner,W. W. (1992) AppliedGroundwaterModeling Simulation oj'Flow and Advective Transport. Academic Press. Christiansen, H., Hill, M. C , Rosbjerg, D. & Jensen, K. H. (1995) Three-dimensional inverse modeling using heads and concentrations at a Danish landfill. In: Models for Assessing Groundwater Quality (ed. by B. J. Wagner, T. H. Illangesekare&K. H. Jensen)(Proc.BoulderSymp.,Colorado,July 1995), 167-175.IAHSPubl.No.227,167175. Hill, M. C. (1992) A computer program (MODFLOWP) for estimating parameters of a transient, three-dimensional groundwater flow model using nonlinear regression. USGS Open File Report 91-484. LeBlanc, D. R. (1984) Digital model of solute transport in a plume of sewage contaminated ground water. In: Movement and fate of solutes in a plume of sewage-contaminatedwater. USGS Open File Report 84-475 (ed. by D. R. LeBlanc) (Cape Cod, Massachusetts). McKenna, S. A. & Poeter, E. P. (1995) Field example of data fusion for site characterization. Wat. Resour. Res. 33(6), 3229-3240. Poeter.E. P. &Hill,M. C. (1996) Inverse models: A necessary next step in groundwater modeling. Groundwater (in press). Poeter, E. P. & McKenna, S. A. (1995) Reducing uncertainty associated with ground-water flow and transport predictions. Groundwater 33(6), 899-904.