- Email: [email protected]
Evaluation of confidence intervals for a steady-state leaky aquifer model
PII: S 0 3 0 9 - 1 7 0 8 ( 9 8 ) 0 0 0 5 5 - 4
Advances in Water Resources Vol. 22, No. 8, pp 807±817, 1999 Ó 1999 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0309-1708/99/$ ± see front matter
Evaluation of con®dence intervals for a steady-state leaky aquifer model S. Christensena,* & R.L. Cooleyb a Department of Earth Sciences, Aarhus University, Ny Munkegade b. 520, 8000 Aarhus C, Denmark Water Resources Division, U.S. Geological Survey, Box 25046, Mail Stop 413, Denver Federal Center, Denver, Colorado 80225, USA
b
(Received 20 May 1997; revised 22 December 1997; accepted 20 November 1998) The fact that dependent variables of groundwater models are generally nonlinear functions of model parameters is shown to be a potentially signi®cant factor in calculating accurate con®dence intervals for both model parameters and functions of the parameters, such as the values of dependent variables calculated by the model. The Lagrangian method of Vecchia and Cooley [Vecchia, A.V. & Cooley, R.L., Water Resources Research, 1987, 23(7), 1237±1250] was used to calculate nonlinear Schee-type con®dence intervals for the parameters and the simulated heads of a steady-state groundwater ¯ow model covering 450 km2 of a leaky aquifer. The nonlinear con®dence intervals are compared to corresponding linear intervals. As suggested by the signi®cant nonlinearity of the regression model, linear con®dence intervals are often not accurate. The commonly made assumption that widths of linear con®dence intervals always underestimate the actual (nonlinear) widths was not correct. Results show that nonlinear eects can cause the nonlinear intervals to be asymmetric and either larger or smaller than the linear approximations. Prior information on transmissivities helps reduce the size of the con®dence intervals, with the most notable eects occurring for the parameters on which there is prior information and for head values in parameter zones for which there is prior information on the parameters. Ó 1999 Elsevier Science Ltd. All rights reserved Key words: con®dence interval, nonlinearity, groundwater ¯ow, model, regression.
model parameters. Cooley and Vecchia13 derived a general method for the calculation of simultaneous con®dence intervals for the output from a groundwater ¯ow model when the statistical distribution of model parameters is known. Vecchia and Cooley23 and Clarke8 independently derived a similar methodology that can be used to compute simultaneous con®dence intervals on both the estimated parameters and the output from a nonlinear regression model with normally distributed residuals. Vecchia and Cooley23 showed that the same algorithm can be used both to estimate the parameters by nonlinear regression and to calculate con®dence limits from which the desired simultaneous con®dence intervals are obtained. Beven and Binley3 presented a Bayesian method that is similar in intent to the method of Vecchia and Cooley23, and Brooks et al.4 presented a nonstatistical variant.
1 INTRODUCTION Estimates of parameters and dependent variables from a calibrated groundwater model are generally uncertain because the data used for calibration are uncertain and because the model never perfectly represents the system or exactly ®ts the data. Con®dence intervals on the estimated (calibrated) model parameters and model dependent variables can be used to express the degree of uncertainty in these quantities, but calculation of con®dence intervals is not straightforward because the solution of the groundwater ¯ow equation for hydraulic heads, and quantities such as ¯ows that are a function of hydraulic heads, is generally a nonlinear function of the *
Corresponding author. 807
808
S. Christensen, R. L. Cooley
Con®dence intervals considered here are Schee-type intervals22, which are simultaneous for all linearizable functions of model parameters and so are the largest of all simultaneous intervals for these functions. Individual con®dence intervals and con®dence intervals that are simultaneous for ®nite numbers of functions, which are smaller than Schee-type intervals, could also have been calculated, but published theory (for example, Ref.16) is not as complete as the theory for the Schee-type intervals (for example, Ref.10). As will be shown, conclusions obtained for the present study will apply to the other types of con®dence intervals. For brevity we will term con®dence intervals calculated using the nonlinear regression model as nonlinear con®dence intervals. The nonlinear regression model can be linearized using a ®rst order Taylor series expansion. Con®dence intervals calculated using the linearized model are the standard ones found in most regression texts22 and are termed here `linear con®dence intervals'. Synthetic case studies9,10,13,17,23 show that: (i) corresponding linear and nonlinear con®dence intervals are often oset or shifted relative to one another, and the nonlinear intervals are generally larger10,17; (ii) the variability in sizes of nonlinear intervals is generally larger than the variability in sizes of corresponding linear intervals10,17; (iii) the dierences between the sizes of corresponding nonlinear and linear con®dence intervals increase as the sizes of the intervals increase9; (iv) use of prior information can have a signi®cant eect on the sizes of con®dence intervals10. In the present study we compute and analyze both linear and nonlinear con®dence intervals using ®eld data to test the calculation procedures and the concepts derived from the synthetic case studies. Linear and nonlinear con®dence intervals are compared for two dierent sets of data used for calibration: (i) observed heads and (ii) observed heads and prior information on the parameters that is obtained from ®eld measurements. 2 BACKGROUND 2.1 Regression model and basis for con®dence intervals The groundwater ¯ow model and accompanying information on the groundwater system are stated as a nonlinear regression model of the form y f b e 1 where y [yi ] is a vector of n observations of the groundwater system (in this study hydraulic head and selected model parameters); b [bj ] is a vector of p unknown, true, model parameters; f(b) [fi (b)] is a vector of n model-computed values corresponding to y; and e [ei ] is a vector of n true errors. The true errors are treated as random variables assumed to have zero means and be distributed normally as
e N 0; Wÿ1 r2 2 where W [Wij ] is an n ´ n known, positive-de®nite weight matrix and r2 is an unknown scalar. Matrix Wÿ1 r2 is the variance±covariance matrix for e. In this study this matrix is assumed to be block diagonal with one block for hydraulic heads and the other for prior information on selected parameters. Diagonal entries of Wÿ1 r2 are variances and o-diagonal entries are covariances. Because all head observations are assumed to be uncorrelated here, all odiagonal entries are zero in the hydraulic head block. The prior information on parameters is assumed to be correlated. Vectors y, f(b), and e are partitioned to conform with the blocks in W. Weights are initially estimated, and can be reestimated as the regression progresses if analysis of residuals indicates that they are incorrect12. Thus, although the weights are assumed to be known, it is apparent that they are approximate. To this extent, the computed con®dence intervals are also approximate. The unknown scalar r2 is estimated by the regression. Schee-type con®dence intervals are to be obtained on some function of parameters g(b). In this study the function g(b) can represent hydraulic head or a parameter, and thus can be the same as f(b). However, in general g(b) can also represent other nonlinear or linear functions of parameters for which a con®dence interval is desired. As indicated earlier, Schee-type intervals are on all linearizable (see Appendix A) functions simultaneously, and thus can be de®ned by the probability statement Prob g bL 6 g b 6 g bU for all linearizable g b 1ÿa
3
where g(bL ) and g(bU ) are the lower and upper con®dence limits, respectively, and sets bL and bU are in general dierent for each function, g(b). These limits are computed as minimum and maximum values of g(b) over the parameter con®dence region. Schee-type intervals should be contrasted with the more common individual con®dence intervals that are on a single function. A good way of visualizing the dierence is to note that, from eqn (3), the probability is a that there will be any (1 ÿ a) ´ 100% Schee-type interval that does not contain the true value, whereas the probability is a that the individual con®dence interval does not contain the true value. Because of this dierence, Schee-type intervals are larger than individual intervals. Calculation of Schee-type con®dence intervals is discussed in Ref.23 and in Appendix A. 2.2 Measures of model nonlinearity Nonlinearity of the regression model results when the model function f(b) is a nonlinear function of model parameters b. There are two components of this non-
Evaluation of con®dence intervals for a steady-state leaky aquifer model linearity, parameter eects nonlinearity and intrinsic nonlinearity2. Parameter eects nonlinearity is the component of the nonlinearity that can in theory be removed by a suitable transformation of model parameters, whereas the intrinsic nonlinearity is the component of the nonlinearity that cannot be removed by any parameter transformation. The sum of the two components is the total nonlinearity, as measured by eqn (4) below, say, and, through some transformation of parameters, it has a minimum achievable value equal to the intrinsic nonlinearity. As described later, nonlinearity aects the Schee-type intervals. Measures of nonlinearity have been developed by Beale2, Linssen20, and Bates and Watts1, among others. Linssen's measures20 are easily computed and are useful to gauge the degree of total and intrinsic nonlinearity of the model. The square of Linssen's measure20 of total nonlinearity can be stated as12 m ÿ P
T ÿ f l ÿ f 0l W f l ÿ f 0l ^ b ps2 l1 M T 2 m P 0 0 ^ ^ fl ÿ f W fl ÿ f
4
l1
where f l f bl ; ^f f ^ b; f 0l is the linear model approximation of f l , computed from ^ bl ÿ ^b 5 f 0l ^f X s2 is the error variance, de®ned as s2 S ^b nÿp
6
and m is the number of parameter sets bl used to com^ is the n ´ p sensitivity pute the measure. In eqn (5) X matrix given by ^ @f i 7 X @bj b^b and in eqn (6) S ^b is the sum of squared errors objective T b. Cooley function for ^b, given by y ÿ f ^ b W y ÿ f ^ 12 and Na give a justi®cation of eqn (4) and recommend using m 2p points bl at the maximum and minimum linear Schee-type con®dence limits for the parameters. From eqn (4) it can be seen that Linssen's measure is a scaled length of the discrepancy between the correct model values fl and the linearized values fol . ^ b for any transformation of The minimum value of M parameters is the intrinsic nonlinearity and is computed using sets of transformed parameters d0l for which the numerator of eqn (4) is a minimum. Linssen20 shows that each set can be found by solving the set of linear equations ^ 0 X ^ T W f l ÿ ^f ^ T WXd X 8 l Thus, the square of Linssen's measure of intrinsic nonlinearity is
m P
W f l ÿ f 0l ^ / ps2 l1 M T 2 m P f 0l ÿ ^f W f 0l ÿ ^f f l ÿ f 0l
809
T
9
l1
where ^ 0 f 0l ^f Xd l
10
Note that fl in eqn (9) is computed using parameter sets bl . Degrees of total nonlinearity have been classi®ed by Beale2 and Cooley and Na12. For example, Beale2 ^b > classi®es a model as disastrously nonlinear if M ^ 1=Fa p; n ÿ p and eectively linear if M b < 0:01= Fa p; n ÿ p, where Fa (p,n ÿ p) is the upper a point of the F distribution with p and n ÿ p degrees of freedom. Cooley and Na12 added an intermediate class ^ b < 0:09=Fa p; n ÿ p, for which linear Schee-type M intervals for parameters were found to approximate the nonlinear ones fairly well in test cases. Large values of ^ b indicate that linear con®dence intervals may not be M accurate. However, small values may not always indi^ b is cate that linear intervals will be accurate because M an average value that may not measure the relevant degree of nonlinearity as it aects a particular interval, and because of the possible in¯uence of nonlinearity of the function g(b). ^ / can be used to correct As shown in Appendix A, M a Schee-type interval for intrinsic nonlinearity. The resulting interval is theoretically conservative to the ^/ degree of approximation used in its derivation, if M 2 were the true value and not an estimate . In practice, the correction should be regarded as approximate and the resulting intervals as probably conservative. 3 FIELD CASE The test case is a steady-state groundwater ¯ow model of a 450 km2 leaky aquifer in Quaternary deposits of glacial till and ¯uvioglacial sand and gravel. The aquifer is located in the western part of the Danish island Zeeland within the catchment of the stream Tude aa (Fig. 1). The aquifer system is outlined in Fig. 2. The aquifer is overlain by 10±80 m of till. In high-lying areas, the aquifer is recharged by downward leakage from a phreatic aquifer in the more permeable upper zone of the till, whereas in low-lying stream valleys and areas near the coast, the aquifer is discharged by upward leakage. More information about the hydrogeology is given by Christensen5. The model was originally calibrated by trial and error5 and subsequently by nonlinear regression7 using MODFLOWP18. Choice of model parameters was based on a detailed description of the hydrogeology5. However, the estimates and linear individual con®dence intervals of several of the 19 groundwater-model
810
S. Christensen, R. L. Cooley
Fig. 3. Ground-water model: (a) zonation and head data (s is estimated standard deviation); (b) simulated head and location of head con®dence intervals.
Fig. 1. Location of aquifer.
Fig. 2. Schematic representation of the aquifer system along the pro®le line in Fig. 1. The arrows give the ¯ow directions in the aquifer and in the overlying till.
parameters in Ref.7 are nearly the same, which indicates redundancy in the parameterization. Deletion of redundant parameters improves the conditioning for the calculation of nonlinear con®dence intervals, so the number of model parameters was reduced to a total of 11: the log10 -transmissivity of the nine zones shown in Fig. 3 Y1 Y9 , and the log10 vertical hydraulic conductivity of the semi-con®ning till layer within zones 1±5 (Z1±5 ) and zones 6±9 (Z6±9 ). Each aquifer zone is characterized by speci®c geological conditions. In zone 1 the aquifer is thin, and well logs show that it is absent in some places. The expected eective transmissivity of this zone is therefore small, which is also indicated by the large head gradients in the northern part of zone 1 (Fig. 3(b)). In the southern part of the zone, the head gradients are smaller because the groundwater ¯uxes are reduced by scattered groundwater withdrawals. The thickness and the transmissivity of the aquifer varies signi®cantly within zones 2, 6 and 7. The aquifer is rather thick and continuous within zones 3±5, with observed transmissivities of the order of
10ÿ3 ±10ÿ2 m2 /s. Zones 8 and 9 represent an extremely heterogeneous area having local sand layers (zone 9) surrounded by till (zone 8). Fig. 3(b) shows that the head gradients are very large in this area. (Note that the observed head ®eld is not shown, but it compares well with the simulated head ®eld in Fig. 3(b).) The simulated leakage is a function of the vertical hydraulic conductivity, the thickness of the semi-con®ning till layer, and of the head dierence between the water table in the upper part of the till and the head in the aquifer. The hydraulic conductivity was estimated by calibration as mentioned, whereas the thickness of the semi-con®ning layer was interpolated from well logs and thus was not a model parameter. The water table in the till was ®xed at an elevation 3 m below the ground level, whereas the head in the aquifer was a dependent variable computed by the model. A no-¯ux condition was used along the entire model boundary. In some areas, the no-¯ux condition is due to geologic boundaries, whereas groundwater divides de®ne the boundary condition in other areas5. Groundwater is withdrawn from several major and minor well ®elds throughout the model area. The total withdrawal is about 8 million m3 per year, and around 40% of the total withdrawal is withdrawn from the four major well ®elds shown in Fig. 1. The withdrawal corresponds to approximately 45% of the ground water leaking to the aquifer, whereas 50% is leaking upward through the till in the low stream valleys, and 5% is leaking to the sea5. 4 PARAMETER ESTIMATION AND MODEL NONLINEARITY The model parameters were estimated by nonlinear regression using a modi®ed code of Cooley and Na12 that uses a combined modi®ed Gauss±Newton and quasiNewton method for the optimization11. The model
Evaluation of con®dence intervals for a steady-state leaky aquifer model parameters were estimated from two dierent sets of data: (a) synchronously measured head in 100 wells (Fig. 3); and (b) these same heads supplemented with prior information on four zonal log10 -transmissivities (Y4 , Y5 , Y6 and Y7 ). The hydraulic heads in wells were observed with an accuracy of a few centimeters. However, due to processes such as small-scale variability in hydrogeology, which are unknown and/or are not represented in the model, the groundwater model can only simulate the large scale variations (here termed ``the drift'') of the true head ®eld7. The weights of the head observations were therefore estimated as the inverses of the variances between a polynomial approximation of the drift6 and the measured heads. The estimated variances vary between 1 and 4.4 m2 (Fig. 3(a)), whereas the covariances are assumed to be zero. Prior information on zonal log-transmissivities was obtained by Christensen7 in the following way. Localscale transmissivity was estimated by analyzing shortduration pumping tests in 90 wells. These results were used to estimate the transmissivity semi-variogram, and to estimate by kriging the local-scale transmissivity at regularly spaced grid points within each of the four zones. The zonal estimates and the corresponding co-
811
variance matrix were obtained by combining the kriged estimates as described by Journel and Huijbregts19 (c.f. pp. 320±324). The inverse covariance matrix was used as the weight matrix of the prior information. The estimated parameters from the nonlinear regression are shown in Fig. 4. The estimated transmissivities lie within the range of transmissivities determined by analysis of numerous pumping tests within the various aquifer zones5,7, and the estimated hydraulic conductivities of the con®ning layer are within the range of values (10ÿ9 ±10ÿ8 m/s) found by analyzing a smaller number of pumping tests5. The parameter estimates also compare well to those of previous model studies in the area5,7. Residual analyses12 show that, in both cases (a) and (b), the weighted residuals, W1=2 e (where e is the residual vector which is an estimate of e), have nearly constant variance; they appear to be unbiased and normally distributed. Thus, there is no reason to suspect that the weight matrix used is not a good estimate of the true covariance matrix for the calibration data (head and prior information). The error variance s2 increases from 1.06 in case (a) to 1.12 in case (b). This may indicate a weak incompatibility between the ¯ow model and the prior information.
Fig. 4. Parameter estimates and con®dence intervals.
812
S. Christensen, R. L. Cooley
^ b is 0.17, and in case In case (a) the total nonlinearity, M ^ (b) M b is 0.12. These values lie between 1/F0:05 (p,n ÿ p) 0.52 and 0.09/F0:05 (p,n ÿ p) 0.05, which means that there is signi®cant nonlinearity in both cases. The large values suggest that linear con®dence intervals may not be ^ / , is 0.06 in case accurate. The intrinsic nonlinearity, M (a) and 0.05 in case (b). Both values are probably not signi®cant, but the con®dence region was corrected to improve the approximation of the con®dence intervals. We used the correction method of Beale2 and Linssen20, which is summarized in Appendix A. Finally note that the prior information slightly reduces the total as well as the intrinsic nonlinearity of the regression model. 5 CALCULATION OF CONFIDENCE INTERVALS In this study, linear and nonlinear 95% Schee-type con®dence intervals were calculated on the estimated parameters and hydraulic head. The linear intervals were calculated by the standard method22, and the nonlinear intervals were calculated by the likelihood method described in Appendix A. The latter intervals were calculated using a modi®ed Gauss±Newton algorithm, as suggested by Vecchia and Cooley23. The calculation of each con®dence limit (the upper or lower bound of the con®dence interval) is a separate optimization problem in which the bound is obtained as an extreme value on the edge of the parameter con®dence region23. To guard against obtaining a local extreme, the calculation of each limit was carried out with nine dierent sets of starting values for parameters. The starting values were obtained by adding or subtracting up to two linearized standard deviations to or from the values estimated by nonlinear regression. This procedure was followed for the calculation of nonlinear intervals for the eleven parameters, as well as for the hydraulic head at 20 grid points homogeneously distributed over the model area (Fig. 3). The variation of parameter starting-values did not produce signi®cant dierences in the con®dence limits for the parameters. Also, for all limits but one, the calculations converged in 5±14 iterations. In the case where no prior information was used, the calculation of the lower bound value of Z1±5 converged more slowly, in 17±25 iterations. The convergence criterion was 0.001 for the relative parameter change from one iteration to the next. The variation of parameter starting values also did not produce signi®cant dierences in the con®dence limits for the head at the grid points. For all limits but one, convergence was achieved in 5±19 iterations. For the upper limit at point 13, however, convergence was not achieved for any of the starting parameters. In this case convergence was achieved by manually adjusting the value of one parameter (Y1 ) and letting the algorithm estimate the other 10 parameter values.
For con®dence intervals on selected parameters and heads, we not only searched for extreme values (i.e. con®dence limits) on the edge of the parameter con®dence region, but also searched the inside of the region. In all of these cases, parameter sets giving the con®dence limits were found on the edge of the region. The results thus strongly indicate that the con®dence intervals are calculated correctly. Based on the above results, we conclude that the algorithm is a fairly robust and ecient method for calculating nonlinear con®dence intervals on both parameters and heads for the present ®eld case. We therefore calculated the head con®dence interval at 485 grid points evenly spaced throughout the model area. At 400 of these points, the calculations converged in less than 25 iterations. These 400 points were used to contour characteristics of the con®dence intervals. The calculation of nonlinear intervals was made on a 50 MHz 80486 PC. The CPU time used was about 5 min per parameter interval, and about 6 min per head interval. 6 RESULTS 6.1 Con®dence intervals on parameters Fig. 4 shows the calculated con®dence intervals for the parameters. For the log10 -vertical conductivities, Z, the sizes of the intervals are about half an order of magnitude, and the linear and nonlinear intervals are quite similar. For the log10 -transmissivities, Y, the sizes of the intervals vary from one to more than three orders of magnitude. The widest intervals are for: (i) the zones with a small Y-estimate and/or small head gradients throughout the zone (Y1 , Y2 , Y4 ), or (ii) a zone without head observations (Y8 ). A comparison of linear and nonlinear con®dence intervals supports conclusions from previous synthetic case studies. Nonlinear intervals are generally larger than the linear intervals, except for Y3 , and the two types of intervals are often oset relative to one another10,17. The dierences among nonlinear intervals are greater than the dierences among corresponding linear intervals10,17. Absolute dierences between the sizes of nonlinear and linear intervals tend to increase as the sizes of intervals increase9, but dierences relative to the sizes of the corresponding linear intervals are more nearly constant. Addition of prior information signi®cantly reduces the sizes of the con®dence intervals for the parameters with prior information, and the linear and nonlinear intervals are similar for these parameters. Intervals for parameters without prior information are nearly unaected by adding prior information on other parameters. Fig. 4 shows that the linear con®dence interval on Y3 is larger than the corresponding nonlinear interval. This contradicts the idea that the linearized variance used to
Evaluation of con®dence intervals for a steady-state leaky aquifer model construct the linear con®dence interval for g(b) corresponds to the Cramer±Rao bound21, and therefore, should be the lower bound. However, the linearized variance only corresponds to the Cramer±Rao bound if it is calculated at the unknown true set of parameters; i.e., for ^b b. Because ^ b in general will dier from b, the linearized variance actually used will generally not correspond to the Cramer±Rao bound, and linear con®dence intervals may therefore be larger than their nonlinear counterparts. 6.2 Con®dence intervals on hydraulic head Fig. 5 shows the calculated con®dence intervals for head at the 20 grid points, and Fig. 6 shows the contoured size of the nonlinear con®dence intervals throughout the model area. The size of the nonlinear intervals generally ranges from 2±10 m within the largest part of the area. The size, however, is as large as 40 m in an area within and east of zones 8 and 9, where head data are lacking. Head gradients in this area are signi®cant (toward two well ®elds in zone 9), which results in a large sensitivity of head to the local parameters (Y7 and Y8 ). Because the
813
estimates of these parameters are uncertain due to the mentioned lack of data, large con®dence intervals result. Some reduction in uncertainty of parameter estimates is expected if head data are collected immediately south and east of zone 8. Similar conditions explain the large con®dence intervals around a minor well ®eld in zone 1. These results show that hydraulic head can vary over a large range when parameter values vary over their con®dence region. Thus, when Schee-type intervals are appropriate, as when gauging the uncertainty of a number of dierent dependent variables and (or) parameters simultaneously, large uncertainties can result. Comparison of the linear and nonlinear intervals in Fig. 5 shows that the two types of intervals are often oset from one another. Many of the nonlinear intervals are skewed, and the skewness tends to increase with the size of the interval. Fig. 7 shows the dierences between the sizes of nonlinear and linear con®dence intervals. Note that the sizes of the nonlinear con®dence intervals in general are larger than the sizes of the linear intervals in the southern part of the model area, whereas the relationship is opposite in the northern part. As argued above,
Fig. 5. Con®dence intervals for the hydraulic head at grid points in Fig. 3(b).
814
S. Christensen, R. L. Cooley
Fig. 6. Size of nonlinear con®dence intervals on head [m].
Fig. 7. Size of nonlinear minus size of linear con®dence interval on head [m].
this indicates that the linearized variance used to compute the linear con®dence interval does not in general correspond to the Cramer±Rao bound. Figs. 5 and 6 show that prior information on Y reduces the width and skewness of the con®dence intervals somewhat in the zones with the prior information, and corresponding linear and nonlinear intervals are more alike in these zones. In contrast to this, the con®dence intervals are nearly unaected in zones without prior information. Fig. 8 shows the reduction in size of the con®dence intervals due to the use of prior information. It can be seen that the prior information reduces the sizes of the intervals by less than 0.5 m in most of the area, but in minor areas the reduction is of the order of 1 m. Within and east of zones 8 and 9, the reduction is, however, up to 16 m. The head gradient in this area is
large, and it is sensitive to Y7 . Prior information signi®cantly reduces the uncertainty of this parameter and thus also reduces the head con®dence interval. 7 SUMMARY AND CONCLUSIONS Nonlinear regression was successfully used to estimate the parameters of a steady-state groundwater ¯ow model of a leaky aquifer. The residuals are normally distributed and nonlinear con®dence intervals can therefore be calculated by the method of Vecchia and Cooley23. Linssen's measure20 shows that total nonlinearity of the regression model is signi®cant, even when prior information on some transmissivities is added to the head data. To a lesser extent, this is also true for the
Evaluation of con®dence intervals for a steady-state leaky aquifer model
815
information. Use of prior information reduced the total nonlinearity of the model only slightly. The present case study compares linear and corresponding nonlinear con®dence intervals and establishes geometric characteristics of nonlinear con®dence intervals. Thus the general conclusions should be valid not only for Schee-type intervals, but also for individual and other types of con®dence intervals as well. ACKNOWLEDGEMENTS This work was partly ®nanced by the Danish Natural Science Research Council via a research fellowship for Steen Christensen. Fig. 8. Reduction of size of nonlinear con®dence interval on head [m] due to the use of prior information on transmissivity.
intrinsic nonlinearity. We therefore corrected the size of the parameter con®dence region used to calculate the nonlinear Schee-type intervals as suggested by Beale2 and Linssen20. The method of Vecchia and Cooley23 was, for the present case, a fairly robust and ecient algorithm for the calculation of nonlinear con®dence intervals on both the hydraulic head and the model parameters. Calculated con®dence limits did not depend on starting parameter values, which indicates that all con®dence limits are probably unique. One of the main assumptions of Vecchia and Cooley's23 Lagrangian method, that the limits of the con®dence interval can be calculated from parameter sets on the edge of the parameter con®dence region, is ful®lled for all the calculated con®dence intervals. As suggested by the signi®cant model nonlinearity, linear con®dence intervals are often not accurate. Nonlinear eects can cause the nonlinear intervals to be oset from, and either larger or smaller than, the linear approximations. The nonlinear head intervals are largest in areas of the model with large head gradients, whereas linear intervals are largest in areas with small gradients. Thus, widths of linear con®dence intervals do not always underestimate the actual (nonlinear) widths so that the linearized variance used to compute a linear con®dence interval does not in general correspond to the Cramer±Rao bound. There is only correspondence if the linearized variance is calculated at the unknown true set of parameters. The general agreement of regression estimates of parameters with and without prior information suggests that the parameter and head data are compatible. The prior information helps reduce and stabilize the con®dence intervals, although the most notable eects occur for the parameters on which there is prior information and for head values in zones for parameters having prior
APPENDIX A It is shown here that if the intrinsic nonlinearity is small, then nearly exact Schee-type con®dence intervals for the function g(b) can be computed by ®nding maximum and minimum values of g(b) over the standard likelihood con®dence region for b. It is also shown that, if the intrinsic nonlinearity is signi®cant, then the size of the con®dence region can be corrected, and the resulting con®dence intervals are conservative. Assume ®rst that intrinsic nonlinearity is negligible. Then there exists some transformation, a, of original model parameters, b, that nearly linearizes the regression model, as measured by eqn (4) 22 (c.f. pp. 133±136). In this case, model function f(b) can be written as f(b) g(a), where g(a) is approximately a linear function of a. In addition the sum of squared errors function S(b) can be written as S(a) using g(a) in place of f(b), and S(b) S(a) 22 (c.f. pp. 195). Because the transformed model is nearly linear, the assumption of normally distributed errors (2) leads to S b ÿ S ^b=p S / ÿ S ^a=p F p; n ÿ p S ^a= n ÿ p S ^b= n ÿ p A:1 where / is the true set of transformed parameters corresponding to the true set b of original parameters, ^a is the set of regression estimates of / corresponding to the set of regression estimates ^b of b, and F(p,n ÿ p) signi®es the F distribution with p and n ÿ p degrees of freedom14 . Based on eqn (A.1), a (1 ÿ a) ´ 100% likelihood con®dence region for b is given by A:2 S b ÿ S ^b 6 ps2 Fa p; n ÿ p da2 where Fa (p,n ÿ p) is the upper a point of F14. As demonstrated subsequently, a Schee-type con®dence interval is given by A:3 min g b; max g b b
b
subject to the constraint
816
S. Christensen, R. L. Cooley
S b ÿ S ^b 6 da2 A:4 If there are no stationary points of g(b) (where og=ob 0) within the con®dence region given by eqn (A.4), or if any stationary point produces an extremum less extreme than one on the boundary of eqn (A.4), then eqn (A.4) can be replaced by the equality constraint S b ÿ S ^ b da2 . In this case the limits given by eqn (A.3) can be computed by ®nding extreme values of h i b A:5 L b; k0 g b k0 da2 ÿ S b S ^ with respect to b and Lagrange multiplier k0 23. An iterative solution for extreme values of eqn (A.5) is given by Vecchia and Cooley23. It is straightforward to show, using a proof analogous to that used by Cooley10 for a Bayesian Scheetype interval, that eqn (A.5) gives a Schee-type con®dence interval for all functions g(b) that are linearizable using Taylor series. However, this can also be demonstrated graphically for p 2 by referring to Fig. 9. To interpret the ®gure, which is a plot in parameter space, note that g(b) constant describes a curve. Also note that g(b1 ), g(b2 ), g(bL ) and g(bU ) are constants. The con®dence region boundary S b ÿ S ^ b da2 is a closed surface (here a closed curve) comprising all possible values for S(b) on the boundary. Now if b actually lies outside of the con®dence region for b, such as b1 , then for at least one function of the form g(b) g(b1 ), such as the one shown, g(b1 ) lies outside of its con®dence interval given by (g(bL ), g(bU )). Thus, the entire curve g(b) g(b1 ) lies outside of the region bounded by g(b) g(bL ) and g(b) g(bU ). However, if b actually lies within the con®dence region, such as b2 , then all linearizable functions of the form g(b) g(b2 ) must lie within their con®dence intervals because they must pass through b2 . Therefore, all linearizable functions are
contained within their con®dence intervals simultaneously with the same relative frequency, 1 ÿ a, that b is contained within its con®dence region. Although it cannot be graphically portrayed, this argument generalizes for any p. We also argue that, because of the physical nature of groundwater models, only linearizable functions are of interest and thus need be considered. If the intrinsic nonlinearity is signi®cant (which can be measured by eqn (9)), then f(b) cannot be nearly linearized and eqn (A.1) may not be nearly correct. As a result, con®dence region eqn (A.2) may not contain b with probability 1 ÿ a. Beale2 gave an approximate correction to obtain a conservative (1 ÿ a) ´ 100% con®dence region. With the correction, da2 is given by n N/ ps2 Fa p; n ÿ p; p 1 da2 1 nÿ1 n p 2 2 N/ ps2 Fa p; n ÿ p; p P 2 da 1 n ÿ pp A:6 where N/ is a theoretical measure of intrinsic nonlinearity obtained by letting sets bl be normally distributed, and then letting m ® 1 in (9)2 . (Actually, Beale2 made an approximation that Guttman and Meeter15 found to yield signi®cant error and Linssen20 corrected. The corrected version is used here.) In this study, we follow ^ / as given by Linssen's20 suggestion and replace N/ by M eqn (9). The signi®cance of the intrinsic nonlinearity is simply measured by the magnitude of the corrected da as compared to the uncorrected value. The same arguments as originally used to justify the use of eqn (A.5) to compute Schee-type intervals can be used to incorporate the correction eqn (A.6). Because the con®dence region is conservative, so are the Scheetype intervals. REFERENCES
Fig. 9. Relations between a con®dence interval (g(bL ), g(bU )) and possible true values g(b1 ) and g(b2 ) for g(b). S b ÿ S ^ b da2 is the con®dence region boundary.
1. Bates, D.M. & Watts, D.G., Relative curvature measures of nonlinearity. Journal of the Royal Statistical Society, Series B, 1980, 42(1), 1±25,. 2. Beale, E.M.L., Con®dence regions in non-linear estimation. Journal of the Royal Statistical Society, series B, 1960, 22(1), 41±76. 3. Beven, K. & Binley, A., The future of distributed models: Model calibration and uncertainty prediction. Hydrological Processes, 1992, 6, 279±298. 4. Brooks, R. J., Lerner, D. N. & Tobias, A. M., Determining the range of predictions of a groundwater model which arises from alternative calibrations. Water Resources Research, 1994, 30(11), 2993±3000. 5. Christensen, S., Hydrological model for the Tude aa catchment. Nordic Hydrology, 1994, 25(3), 145±166. 6. Christensen, S., Calibration data accuracy in a regionalscale groundwater model. Internal research report, Dept. of Earth Sciences, Aarhus University, Aarhus, 1995, 29 pp. 7. Christensen, S., On the strategy of estimating regionalscale transmissivity ®elds. Ground Water, 1997, 35(1), 131± 139.
Evaluation of con®dence intervals for a steady-state leaky aquifer model 8. Clarke, G.P.Y., Approximate con®dence limits for a parameter function in nonlinear regression. Journal of the American Statistical Association, 1987, 82(397), 221±230. 9. Cooley, R.L., Exact Schee-type con®dence intervals for output from groundwater ¯ow models, 1, Use of hydrogeologic information. Water Resources Research, 1993, 29(1), 17±33. 10. Cooley, R.L., Exact Schee-type con®dence intervals for output from groundwater ¯ow models, 2, Combined use of hydrogeologic information and calibration data. Water Resources Research, 1993, 29(1), 35±50. 11. Cooley, R.L. & Hill, M.C., A comparison of three Newton-like nonlinear least-squares methods for estimating parameters of ground-water ¯ow models. Computational Methods in Water Resources IX, Vol. 1, (ed. by T.F. Russell, R.E. Ewing, C.A. Brebbia, W.G. Gray & G.F. Pinder), 1992, 379±386. 12. Cooley, R.L. & Na, R.L., Regression modelling of ground-water ¯ow. U.S. Geological Survey Techniques of Water-Resources Investigations, 1990, book 3, ch. B4, 232 pp. 13. Cooley, R.L. & Vecchia, A.V., Calculation of nonlinear con®dence and prediction intervals for ground-water ¯ow models. Water Resources Bulletin, 1987, 23(4), 581±599. 14. Draper, N. & Smith, H., Applied regression analysis. 2nd ed., John Wiley, New York, 1981, 709 pp.
817
15. Guttman, I. & Meeter, D.A., On Beale's measures of nonlinearity. Technometrics, 1965, 7(4), 623±637. 16. Hamilton, D., & Wiens, D., Correction factors for F ratios in nonlinear regression. Biometrica, 1987, 74(2), 423±425. 17. Hill, M.C., Analysis of accuracy of approximate, simultaneous, nonlinear con®dence intervals on hydraulic heads in analytical and numerical test cases. Water Resources Research, 1989, 25(21), 177±190. 18. Hill, M.C., A computer program (MODFLOWP) for estimating parameters of a transient, three-dimensional ground-water ¯ow model using nonlinear regression. U.S. Geological Survey, Open-File Report 91-484, 1992, 358 pp. 19. Journel, A.G. & Huijbregts, C., Mining geostatistics. Academic Press, New York, 1978, 600 pp. 20. Linssen, H.N., Nonlinearity measures: A case study. Statistica Neerlandica, 1975, 29, 93±99. 21. Mood, A.M., Graybill, F.A. & Boes, D.C., Introduction to the theory of statistics. 3rd ed., McGraw-Hill, New York, 1974, 564 pp. 22. Seber, G.A.F. & Wild, C.J., Nonlinear regression. John Wiley, New York, 1989, 768 pp. 23. Vecchia, A.V. & Cooley, R.L., Simultaneous con®dence and prediction intervals for nonlinear regression models with application to a groundwater ¯ow model. Water Resources Research, 1987, 23(7), 1237±1250.
Advances in Water Resources Vol. 22, No. 8, pp 807±817, 1999 Ó 1999 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0309-1708/99/$ ± see front matter
Evaluation of con®dence intervals for a steady-state leaky aquifer model S. Christensena,* & R.L. Cooleyb a Department of Earth Sciences, Aarhus University, Ny Munkegade b. 520, 8000 Aarhus C, Denmark Water Resources Division, U.S. Geological Survey, Box 25046, Mail Stop 413, Denver Federal Center, Denver, Colorado 80225, USA
b
(Received 20 May 1997; revised 22 December 1997; accepted 20 November 1998) The fact that dependent variables of groundwater models are generally nonlinear functions of model parameters is shown to be a potentially signi®cant factor in calculating accurate con®dence intervals for both model parameters and functions of the parameters, such as the values of dependent variables calculated by the model. The Lagrangian method of Vecchia and Cooley [Vecchia, A.V. & Cooley, R.L., Water Resources Research, 1987, 23(7), 1237±1250] was used to calculate nonlinear Schee-type con®dence intervals for the parameters and the simulated heads of a steady-state groundwater ¯ow model covering 450 km2 of a leaky aquifer. The nonlinear con®dence intervals are compared to corresponding linear intervals. As suggested by the signi®cant nonlinearity of the regression model, linear con®dence intervals are often not accurate. The commonly made assumption that widths of linear con®dence intervals always underestimate the actual (nonlinear) widths was not correct. Results show that nonlinear eects can cause the nonlinear intervals to be asymmetric and either larger or smaller than the linear approximations. Prior information on transmissivities helps reduce the size of the con®dence intervals, with the most notable eects occurring for the parameters on which there is prior information and for head values in parameter zones for which there is prior information on the parameters. Ó 1999 Elsevier Science Ltd. All rights reserved Key words: con®dence interval, nonlinearity, groundwater ¯ow, model, regression.
model parameters. Cooley and Vecchia13 derived a general method for the calculation of simultaneous con®dence intervals for the output from a groundwater ¯ow model when the statistical distribution of model parameters is known. Vecchia and Cooley23 and Clarke8 independently derived a similar methodology that can be used to compute simultaneous con®dence intervals on both the estimated parameters and the output from a nonlinear regression model with normally distributed residuals. Vecchia and Cooley23 showed that the same algorithm can be used both to estimate the parameters by nonlinear regression and to calculate con®dence limits from which the desired simultaneous con®dence intervals are obtained. Beven and Binley3 presented a Bayesian method that is similar in intent to the method of Vecchia and Cooley23, and Brooks et al.4 presented a nonstatistical variant.
1 INTRODUCTION Estimates of parameters and dependent variables from a calibrated groundwater model are generally uncertain because the data used for calibration are uncertain and because the model never perfectly represents the system or exactly ®ts the data. Con®dence intervals on the estimated (calibrated) model parameters and model dependent variables can be used to express the degree of uncertainty in these quantities, but calculation of con®dence intervals is not straightforward because the solution of the groundwater ¯ow equation for hydraulic heads, and quantities such as ¯ows that are a function of hydraulic heads, is generally a nonlinear function of the *
Corresponding author. 807
808
S. Christensen, R. L. Cooley
Con®dence intervals considered here are Schee-type intervals22, which are simultaneous for all linearizable functions of model parameters and so are the largest of all simultaneous intervals for these functions. Individual con®dence intervals and con®dence intervals that are simultaneous for ®nite numbers of functions, which are smaller than Schee-type intervals, could also have been calculated, but published theory (for example, Ref.16) is not as complete as the theory for the Schee-type intervals (for example, Ref.10). As will be shown, conclusions obtained for the present study will apply to the other types of con®dence intervals. For brevity we will term con®dence intervals calculated using the nonlinear regression model as nonlinear con®dence intervals. The nonlinear regression model can be linearized using a ®rst order Taylor series expansion. Con®dence intervals calculated using the linearized model are the standard ones found in most regression texts22 and are termed here `linear con®dence intervals'. Synthetic case studies9,10,13,17,23 show that: (i) corresponding linear and nonlinear con®dence intervals are often oset or shifted relative to one another, and the nonlinear intervals are generally larger10,17; (ii) the variability in sizes of nonlinear intervals is generally larger than the variability in sizes of corresponding linear intervals10,17; (iii) the dierences between the sizes of corresponding nonlinear and linear con®dence intervals increase as the sizes of the intervals increase9; (iv) use of prior information can have a signi®cant eect on the sizes of con®dence intervals10. In the present study we compute and analyze both linear and nonlinear con®dence intervals using ®eld data to test the calculation procedures and the concepts derived from the synthetic case studies. Linear and nonlinear con®dence intervals are compared for two dierent sets of data used for calibration: (i) observed heads and (ii) observed heads and prior information on the parameters that is obtained from ®eld measurements. 2 BACKGROUND 2.1 Regression model and basis for con®dence intervals The groundwater ¯ow model and accompanying information on the groundwater system are stated as a nonlinear regression model of the form y f b e 1 where y [yi ] is a vector of n observations of the groundwater system (in this study hydraulic head and selected model parameters); b [bj ] is a vector of p unknown, true, model parameters; f(b) [fi (b)] is a vector of n model-computed values corresponding to y; and e [ei ] is a vector of n true errors. The true errors are treated as random variables assumed to have zero means and be distributed normally as
e N 0; Wÿ1 r2 2 where W [Wij ] is an n ´ n known, positive-de®nite weight matrix and r2 is an unknown scalar. Matrix Wÿ1 r2 is the variance±covariance matrix for e. In this study this matrix is assumed to be block diagonal with one block for hydraulic heads and the other for prior information on selected parameters. Diagonal entries of Wÿ1 r2 are variances and o-diagonal entries are covariances. Because all head observations are assumed to be uncorrelated here, all odiagonal entries are zero in the hydraulic head block. The prior information on parameters is assumed to be correlated. Vectors y, f(b), and e are partitioned to conform with the blocks in W. Weights are initially estimated, and can be reestimated as the regression progresses if analysis of residuals indicates that they are incorrect12. Thus, although the weights are assumed to be known, it is apparent that they are approximate. To this extent, the computed con®dence intervals are also approximate. The unknown scalar r2 is estimated by the regression. Schee-type con®dence intervals are to be obtained on some function of parameters g(b). In this study the function g(b) can represent hydraulic head or a parameter, and thus can be the same as f(b). However, in general g(b) can also represent other nonlinear or linear functions of parameters for which a con®dence interval is desired. As indicated earlier, Schee-type intervals are on all linearizable (see Appendix A) functions simultaneously, and thus can be de®ned by the probability statement Prob g bL 6 g b 6 g bU for all linearizable g b 1ÿa
3
where g(bL ) and g(bU ) are the lower and upper con®dence limits, respectively, and sets bL and bU are in general dierent for each function, g(b). These limits are computed as minimum and maximum values of g(b) over the parameter con®dence region. Schee-type intervals should be contrasted with the more common individual con®dence intervals that are on a single function. A good way of visualizing the dierence is to note that, from eqn (3), the probability is a that there will be any (1 ÿ a) ´ 100% Schee-type interval that does not contain the true value, whereas the probability is a that the individual con®dence interval does not contain the true value. Because of this dierence, Schee-type intervals are larger than individual intervals. Calculation of Schee-type con®dence intervals is discussed in Ref.23 and in Appendix A. 2.2 Measures of model nonlinearity Nonlinearity of the regression model results when the model function f(b) is a nonlinear function of model parameters b. There are two components of this non-
Evaluation of con®dence intervals for a steady-state leaky aquifer model linearity, parameter eects nonlinearity and intrinsic nonlinearity2. Parameter eects nonlinearity is the component of the nonlinearity that can in theory be removed by a suitable transformation of model parameters, whereas the intrinsic nonlinearity is the component of the nonlinearity that cannot be removed by any parameter transformation. The sum of the two components is the total nonlinearity, as measured by eqn (4) below, say, and, through some transformation of parameters, it has a minimum achievable value equal to the intrinsic nonlinearity. As described later, nonlinearity aects the Schee-type intervals. Measures of nonlinearity have been developed by Beale2, Linssen20, and Bates and Watts1, among others. Linssen's measures20 are easily computed and are useful to gauge the degree of total and intrinsic nonlinearity of the model. The square of Linssen's measure20 of total nonlinearity can be stated as12 m ÿ P
T ÿ f l ÿ f 0l W f l ÿ f 0l ^ b ps2 l1 M T 2 m P 0 0 ^ ^ fl ÿ f W fl ÿ f
4
l1
where f l f bl ; ^f f ^ b; f 0l is the linear model approximation of f l , computed from ^ bl ÿ ^b 5 f 0l ^f X s2 is the error variance, de®ned as s2 S ^b nÿp
6
and m is the number of parameter sets bl used to com^ is the n ´ p sensitivity pute the measure. In eqn (5) X matrix given by ^ @f i 7 X @bj b^b and in eqn (6) S ^b is the sum of squared errors objective T b. Cooley function for ^b, given by y ÿ f ^ b W y ÿ f ^ 12 and Na give a justi®cation of eqn (4) and recommend using m 2p points bl at the maximum and minimum linear Schee-type con®dence limits for the parameters. From eqn (4) it can be seen that Linssen's measure is a scaled length of the discrepancy between the correct model values fl and the linearized values fol . ^ b for any transformation of The minimum value of M parameters is the intrinsic nonlinearity and is computed using sets of transformed parameters d0l for which the numerator of eqn (4) is a minimum. Linssen20 shows that each set can be found by solving the set of linear equations ^ 0 X ^ T W f l ÿ ^f ^ T WXd X 8 l Thus, the square of Linssen's measure of intrinsic nonlinearity is
m P
W f l ÿ f 0l ^ / ps2 l1 M T 2 m P f 0l ÿ ^f W f 0l ÿ ^f f l ÿ f 0l
809
T
9
l1
where ^ 0 f 0l ^f Xd l
10
Note that fl in eqn (9) is computed using parameter sets bl . Degrees of total nonlinearity have been classi®ed by Beale2 and Cooley and Na12. For example, Beale2 ^b > classi®es a model as disastrously nonlinear if M ^ 1=Fa p; n ÿ p and eectively linear if M b < 0:01= Fa p; n ÿ p, where Fa (p,n ÿ p) is the upper a point of the F distribution with p and n ÿ p degrees of freedom. Cooley and Na12 added an intermediate class ^ b < 0:09=Fa p; n ÿ p, for which linear Schee-type M intervals for parameters were found to approximate the nonlinear ones fairly well in test cases. Large values of ^ b indicate that linear con®dence intervals may not be M accurate. However, small values may not always indi^ b is cate that linear intervals will be accurate because M an average value that may not measure the relevant degree of nonlinearity as it aects a particular interval, and because of the possible in¯uence of nonlinearity of the function g(b). ^ / can be used to correct As shown in Appendix A, M a Schee-type interval for intrinsic nonlinearity. The resulting interval is theoretically conservative to the ^/ degree of approximation used in its derivation, if M 2 were the true value and not an estimate . In practice, the correction should be regarded as approximate and the resulting intervals as probably conservative. 3 FIELD CASE The test case is a steady-state groundwater ¯ow model of a 450 km2 leaky aquifer in Quaternary deposits of glacial till and ¯uvioglacial sand and gravel. The aquifer is located in the western part of the Danish island Zeeland within the catchment of the stream Tude aa (Fig. 1). The aquifer system is outlined in Fig. 2. The aquifer is overlain by 10±80 m of till. In high-lying areas, the aquifer is recharged by downward leakage from a phreatic aquifer in the more permeable upper zone of the till, whereas in low-lying stream valleys and areas near the coast, the aquifer is discharged by upward leakage. More information about the hydrogeology is given by Christensen5. The model was originally calibrated by trial and error5 and subsequently by nonlinear regression7 using MODFLOWP18. Choice of model parameters was based on a detailed description of the hydrogeology5. However, the estimates and linear individual con®dence intervals of several of the 19 groundwater-model
810
S. Christensen, R. L. Cooley
Fig. 3. Ground-water model: (a) zonation and head data (s is estimated standard deviation); (b) simulated head and location of head con®dence intervals.
Fig. 1. Location of aquifer.
Fig. 2. Schematic representation of the aquifer system along the pro®le line in Fig. 1. The arrows give the ¯ow directions in the aquifer and in the overlying till.
parameters in Ref.7 are nearly the same, which indicates redundancy in the parameterization. Deletion of redundant parameters improves the conditioning for the calculation of nonlinear con®dence intervals, so the number of model parameters was reduced to a total of 11: the log10 -transmissivity of the nine zones shown in Fig. 3 Y1 Y9 , and the log10 vertical hydraulic conductivity of the semi-con®ning till layer within zones 1±5 (Z1±5 ) and zones 6±9 (Z6±9 ). Each aquifer zone is characterized by speci®c geological conditions. In zone 1 the aquifer is thin, and well logs show that it is absent in some places. The expected eective transmissivity of this zone is therefore small, which is also indicated by the large head gradients in the northern part of zone 1 (Fig. 3(b)). In the southern part of the zone, the head gradients are smaller because the groundwater ¯uxes are reduced by scattered groundwater withdrawals. The thickness and the transmissivity of the aquifer varies signi®cantly within zones 2, 6 and 7. The aquifer is rather thick and continuous within zones 3±5, with observed transmissivities of the order of
10ÿ3 ±10ÿ2 m2 /s. Zones 8 and 9 represent an extremely heterogeneous area having local sand layers (zone 9) surrounded by till (zone 8). Fig. 3(b) shows that the head gradients are very large in this area. (Note that the observed head ®eld is not shown, but it compares well with the simulated head ®eld in Fig. 3(b).) The simulated leakage is a function of the vertical hydraulic conductivity, the thickness of the semi-con®ning till layer, and of the head dierence between the water table in the upper part of the till and the head in the aquifer. The hydraulic conductivity was estimated by calibration as mentioned, whereas the thickness of the semi-con®ning layer was interpolated from well logs and thus was not a model parameter. The water table in the till was ®xed at an elevation 3 m below the ground level, whereas the head in the aquifer was a dependent variable computed by the model. A no-¯ux condition was used along the entire model boundary. In some areas, the no-¯ux condition is due to geologic boundaries, whereas groundwater divides de®ne the boundary condition in other areas5. Groundwater is withdrawn from several major and minor well ®elds throughout the model area. The total withdrawal is about 8 million m3 per year, and around 40% of the total withdrawal is withdrawn from the four major well ®elds shown in Fig. 1. The withdrawal corresponds to approximately 45% of the ground water leaking to the aquifer, whereas 50% is leaking upward through the till in the low stream valleys, and 5% is leaking to the sea5. 4 PARAMETER ESTIMATION AND MODEL NONLINEARITY The model parameters were estimated by nonlinear regression using a modi®ed code of Cooley and Na12 that uses a combined modi®ed Gauss±Newton and quasiNewton method for the optimization11. The model
Evaluation of con®dence intervals for a steady-state leaky aquifer model parameters were estimated from two dierent sets of data: (a) synchronously measured head in 100 wells (Fig. 3); and (b) these same heads supplemented with prior information on four zonal log10 -transmissivities (Y4 , Y5 , Y6 and Y7 ). The hydraulic heads in wells were observed with an accuracy of a few centimeters. However, due to processes such as small-scale variability in hydrogeology, which are unknown and/or are not represented in the model, the groundwater model can only simulate the large scale variations (here termed ``the drift'') of the true head ®eld7. The weights of the head observations were therefore estimated as the inverses of the variances between a polynomial approximation of the drift6 and the measured heads. The estimated variances vary between 1 and 4.4 m2 (Fig. 3(a)), whereas the covariances are assumed to be zero. Prior information on zonal log-transmissivities was obtained by Christensen7 in the following way. Localscale transmissivity was estimated by analyzing shortduration pumping tests in 90 wells. These results were used to estimate the transmissivity semi-variogram, and to estimate by kriging the local-scale transmissivity at regularly spaced grid points within each of the four zones. The zonal estimates and the corresponding co-
811
variance matrix were obtained by combining the kriged estimates as described by Journel and Huijbregts19 (c.f. pp. 320±324). The inverse covariance matrix was used as the weight matrix of the prior information. The estimated parameters from the nonlinear regression are shown in Fig. 4. The estimated transmissivities lie within the range of transmissivities determined by analysis of numerous pumping tests within the various aquifer zones5,7, and the estimated hydraulic conductivities of the con®ning layer are within the range of values (10ÿ9 ±10ÿ8 m/s) found by analyzing a smaller number of pumping tests5. The parameter estimates also compare well to those of previous model studies in the area5,7. Residual analyses12 show that, in both cases (a) and (b), the weighted residuals, W1=2 e (where e is the residual vector which is an estimate of e), have nearly constant variance; they appear to be unbiased and normally distributed. Thus, there is no reason to suspect that the weight matrix used is not a good estimate of the true covariance matrix for the calibration data (head and prior information). The error variance s2 increases from 1.06 in case (a) to 1.12 in case (b). This may indicate a weak incompatibility between the ¯ow model and the prior information.
Fig. 4. Parameter estimates and con®dence intervals.
812
S. Christensen, R. L. Cooley
^ b is 0.17, and in case In case (a) the total nonlinearity, M ^ (b) M b is 0.12. These values lie between 1/F0:05 (p,n ÿ p) 0.52 and 0.09/F0:05 (p,n ÿ p) 0.05, which means that there is signi®cant nonlinearity in both cases. The large values suggest that linear con®dence intervals may not be ^ / , is 0.06 in case accurate. The intrinsic nonlinearity, M (a) and 0.05 in case (b). Both values are probably not signi®cant, but the con®dence region was corrected to improve the approximation of the con®dence intervals. We used the correction method of Beale2 and Linssen20, which is summarized in Appendix A. Finally note that the prior information slightly reduces the total as well as the intrinsic nonlinearity of the regression model. 5 CALCULATION OF CONFIDENCE INTERVALS In this study, linear and nonlinear 95% Schee-type con®dence intervals were calculated on the estimated parameters and hydraulic head. The linear intervals were calculated by the standard method22, and the nonlinear intervals were calculated by the likelihood method described in Appendix A. The latter intervals were calculated using a modi®ed Gauss±Newton algorithm, as suggested by Vecchia and Cooley23. The calculation of each con®dence limit (the upper or lower bound of the con®dence interval) is a separate optimization problem in which the bound is obtained as an extreme value on the edge of the parameter con®dence region23. To guard against obtaining a local extreme, the calculation of each limit was carried out with nine dierent sets of starting values for parameters. The starting values were obtained by adding or subtracting up to two linearized standard deviations to or from the values estimated by nonlinear regression. This procedure was followed for the calculation of nonlinear intervals for the eleven parameters, as well as for the hydraulic head at 20 grid points homogeneously distributed over the model area (Fig. 3). The variation of parameter starting-values did not produce signi®cant dierences in the con®dence limits for the parameters. Also, for all limits but one, the calculations converged in 5±14 iterations. In the case where no prior information was used, the calculation of the lower bound value of Z1±5 converged more slowly, in 17±25 iterations. The convergence criterion was 0.001 for the relative parameter change from one iteration to the next. The variation of parameter starting values also did not produce signi®cant dierences in the con®dence limits for the head at the grid points. For all limits but one, convergence was achieved in 5±19 iterations. For the upper limit at point 13, however, convergence was not achieved for any of the starting parameters. In this case convergence was achieved by manually adjusting the value of one parameter (Y1 ) and letting the algorithm estimate the other 10 parameter values.
For con®dence intervals on selected parameters and heads, we not only searched for extreme values (i.e. con®dence limits) on the edge of the parameter con®dence region, but also searched the inside of the region. In all of these cases, parameter sets giving the con®dence limits were found on the edge of the region. The results thus strongly indicate that the con®dence intervals are calculated correctly. Based on the above results, we conclude that the algorithm is a fairly robust and ecient method for calculating nonlinear con®dence intervals on both parameters and heads for the present ®eld case. We therefore calculated the head con®dence interval at 485 grid points evenly spaced throughout the model area. At 400 of these points, the calculations converged in less than 25 iterations. These 400 points were used to contour characteristics of the con®dence intervals. The calculation of nonlinear intervals was made on a 50 MHz 80486 PC. The CPU time used was about 5 min per parameter interval, and about 6 min per head interval. 6 RESULTS 6.1 Con®dence intervals on parameters Fig. 4 shows the calculated con®dence intervals for the parameters. For the log10 -vertical conductivities, Z, the sizes of the intervals are about half an order of magnitude, and the linear and nonlinear intervals are quite similar. For the log10 -transmissivities, Y, the sizes of the intervals vary from one to more than three orders of magnitude. The widest intervals are for: (i) the zones with a small Y-estimate and/or small head gradients throughout the zone (Y1 , Y2 , Y4 ), or (ii) a zone without head observations (Y8 ). A comparison of linear and nonlinear con®dence intervals supports conclusions from previous synthetic case studies. Nonlinear intervals are generally larger than the linear intervals, except for Y3 , and the two types of intervals are often oset relative to one another10,17. The dierences among nonlinear intervals are greater than the dierences among corresponding linear intervals10,17. Absolute dierences between the sizes of nonlinear and linear intervals tend to increase as the sizes of intervals increase9, but dierences relative to the sizes of the corresponding linear intervals are more nearly constant. Addition of prior information signi®cantly reduces the sizes of the con®dence intervals for the parameters with prior information, and the linear and nonlinear intervals are similar for these parameters. Intervals for parameters without prior information are nearly unaected by adding prior information on other parameters. Fig. 4 shows that the linear con®dence interval on Y3 is larger than the corresponding nonlinear interval. This contradicts the idea that the linearized variance used to
Evaluation of con®dence intervals for a steady-state leaky aquifer model construct the linear con®dence interval for g(b) corresponds to the Cramer±Rao bound21, and therefore, should be the lower bound. However, the linearized variance only corresponds to the Cramer±Rao bound if it is calculated at the unknown true set of parameters; i.e., for ^b b. Because ^ b in general will dier from b, the linearized variance actually used will generally not correspond to the Cramer±Rao bound, and linear con®dence intervals may therefore be larger than their nonlinear counterparts. 6.2 Con®dence intervals on hydraulic head Fig. 5 shows the calculated con®dence intervals for head at the 20 grid points, and Fig. 6 shows the contoured size of the nonlinear con®dence intervals throughout the model area. The size of the nonlinear intervals generally ranges from 2±10 m within the largest part of the area. The size, however, is as large as 40 m in an area within and east of zones 8 and 9, where head data are lacking. Head gradients in this area are signi®cant (toward two well ®elds in zone 9), which results in a large sensitivity of head to the local parameters (Y7 and Y8 ). Because the
813
estimates of these parameters are uncertain due to the mentioned lack of data, large con®dence intervals result. Some reduction in uncertainty of parameter estimates is expected if head data are collected immediately south and east of zone 8. Similar conditions explain the large con®dence intervals around a minor well ®eld in zone 1. These results show that hydraulic head can vary over a large range when parameter values vary over their con®dence region. Thus, when Schee-type intervals are appropriate, as when gauging the uncertainty of a number of dierent dependent variables and (or) parameters simultaneously, large uncertainties can result. Comparison of the linear and nonlinear intervals in Fig. 5 shows that the two types of intervals are often oset from one another. Many of the nonlinear intervals are skewed, and the skewness tends to increase with the size of the interval. Fig. 7 shows the dierences between the sizes of nonlinear and linear con®dence intervals. Note that the sizes of the nonlinear con®dence intervals in general are larger than the sizes of the linear intervals in the southern part of the model area, whereas the relationship is opposite in the northern part. As argued above,
Fig. 5. Con®dence intervals for the hydraulic head at grid points in Fig. 3(b).
814
S. Christensen, R. L. Cooley
Fig. 6. Size of nonlinear con®dence intervals on head [m].
Fig. 7. Size of nonlinear minus size of linear con®dence interval on head [m].
this indicates that the linearized variance used to compute the linear con®dence interval does not in general correspond to the Cramer±Rao bound. Figs. 5 and 6 show that prior information on Y reduces the width and skewness of the con®dence intervals somewhat in the zones with the prior information, and corresponding linear and nonlinear intervals are more alike in these zones. In contrast to this, the con®dence intervals are nearly unaected in zones without prior information. Fig. 8 shows the reduction in size of the con®dence intervals due to the use of prior information. It can be seen that the prior information reduces the sizes of the intervals by less than 0.5 m in most of the area, but in minor areas the reduction is of the order of 1 m. Within and east of zones 8 and 9, the reduction is, however, up to 16 m. The head gradient in this area is
large, and it is sensitive to Y7 . Prior information signi®cantly reduces the uncertainty of this parameter and thus also reduces the head con®dence interval. 7 SUMMARY AND CONCLUSIONS Nonlinear regression was successfully used to estimate the parameters of a steady-state groundwater ¯ow model of a leaky aquifer. The residuals are normally distributed and nonlinear con®dence intervals can therefore be calculated by the method of Vecchia and Cooley23. Linssen's measure20 shows that total nonlinearity of the regression model is signi®cant, even when prior information on some transmissivities is added to the head data. To a lesser extent, this is also true for the
Evaluation of con®dence intervals for a steady-state leaky aquifer model
815
information. Use of prior information reduced the total nonlinearity of the model only slightly. The present case study compares linear and corresponding nonlinear con®dence intervals and establishes geometric characteristics of nonlinear con®dence intervals. Thus the general conclusions should be valid not only for Schee-type intervals, but also for individual and other types of con®dence intervals as well. ACKNOWLEDGEMENTS This work was partly ®nanced by the Danish Natural Science Research Council via a research fellowship for Steen Christensen. Fig. 8. Reduction of size of nonlinear con®dence interval on head [m] due to the use of prior information on transmissivity.
intrinsic nonlinearity. We therefore corrected the size of the parameter con®dence region used to calculate the nonlinear Schee-type intervals as suggested by Beale2 and Linssen20. The method of Vecchia and Cooley23 was, for the present case, a fairly robust and ecient algorithm for the calculation of nonlinear con®dence intervals on both the hydraulic head and the model parameters. Calculated con®dence limits did not depend on starting parameter values, which indicates that all con®dence limits are probably unique. One of the main assumptions of Vecchia and Cooley's23 Lagrangian method, that the limits of the con®dence interval can be calculated from parameter sets on the edge of the parameter con®dence region, is ful®lled for all the calculated con®dence intervals. As suggested by the signi®cant model nonlinearity, linear con®dence intervals are often not accurate. Nonlinear eects can cause the nonlinear intervals to be oset from, and either larger or smaller than, the linear approximations. The nonlinear head intervals are largest in areas of the model with large head gradients, whereas linear intervals are largest in areas with small gradients. Thus, widths of linear con®dence intervals do not always underestimate the actual (nonlinear) widths so that the linearized variance used to compute a linear con®dence interval does not in general correspond to the Cramer±Rao bound. There is only correspondence if the linearized variance is calculated at the unknown true set of parameters. The general agreement of regression estimates of parameters with and without prior information suggests that the parameter and head data are compatible. The prior information helps reduce and stabilize the con®dence intervals, although the most notable eects occur for the parameters on which there is prior information and for head values in zones for parameters having prior
APPENDIX A It is shown here that if the intrinsic nonlinearity is small, then nearly exact Schee-type con®dence intervals for the function g(b) can be computed by ®nding maximum and minimum values of g(b) over the standard likelihood con®dence region for b. It is also shown that, if the intrinsic nonlinearity is signi®cant, then the size of the con®dence region can be corrected, and the resulting con®dence intervals are conservative. Assume ®rst that intrinsic nonlinearity is negligible. Then there exists some transformation, a, of original model parameters, b, that nearly linearizes the regression model, as measured by eqn (4) 22 (c.f. pp. 133±136). In this case, model function f(b) can be written as f(b) g(a), where g(a) is approximately a linear function of a. In addition the sum of squared errors function S(b) can be written as S(a) using g(a) in place of f(b), and S(b) S(a) 22 (c.f. pp. 195). Because the transformed model is nearly linear, the assumption of normally distributed errors (2) leads to S b ÿ S ^b=p S / ÿ S ^a=p F p; n ÿ p S ^a= n ÿ p S ^b= n ÿ p A:1 where / is the true set of transformed parameters corresponding to the true set b of original parameters, ^a is the set of regression estimates of / corresponding to the set of regression estimates ^b of b, and F(p,n ÿ p) signi®es the F distribution with p and n ÿ p degrees of freedom14 . Based on eqn (A.1), a (1 ÿ a) ´ 100% likelihood con®dence region for b is given by A:2 S b ÿ S ^b 6 ps2 Fa p; n ÿ p da2 where Fa (p,n ÿ p) is the upper a point of F14. As demonstrated subsequently, a Schee-type con®dence interval is given by A:3 min g b; max g b b
b
subject to the constraint
816
S. Christensen, R. L. Cooley
S b ÿ S ^b 6 da2 A:4 If there are no stationary points of g(b) (where og=ob 0) within the con®dence region given by eqn (A.4), or if any stationary point produces an extremum less extreme than one on the boundary of eqn (A.4), then eqn (A.4) can be replaced by the equality constraint S b ÿ S ^ b da2 . In this case the limits given by eqn (A.3) can be computed by ®nding extreme values of h i b A:5 L b; k0 g b k0 da2 ÿ S b S ^ with respect to b and Lagrange multiplier k0 23. An iterative solution for extreme values of eqn (A.5) is given by Vecchia and Cooley23. It is straightforward to show, using a proof analogous to that used by Cooley10 for a Bayesian Scheetype interval, that eqn (A.5) gives a Schee-type con®dence interval for all functions g(b) that are linearizable using Taylor series. However, this can also be demonstrated graphically for p 2 by referring to Fig. 9. To interpret the ®gure, which is a plot in parameter space, note that g(b) constant describes a curve. Also note that g(b1 ), g(b2 ), g(bL ) and g(bU ) are constants. The con®dence region boundary S b ÿ S ^ b da2 is a closed surface (here a closed curve) comprising all possible values for S(b) on the boundary. Now if b actually lies outside of the con®dence region for b, such as b1 , then for at least one function of the form g(b) g(b1 ), such as the one shown, g(b1 ) lies outside of its con®dence interval given by (g(bL ), g(bU )). Thus, the entire curve g(b) g(b1 ) lies outside of the region bounded by g(b) g(bL ) and g(b) g(bU ). However, if b actually lies within the con®dence region, such as b2 , then all linearizable functions of the form g(b) g(b2 ) must lie within their con®dence intervals because they must pass through b2 . Therefore, all linearizable functions are
contained within their con®dence intervals simultaneously with the same relative frequency, 1 ÿ a, that b is contained within its con®dence region. Although it cannot be graphically portrayed, this argument generalizes for any p. We also argue that, because of the physical nature of groundwater models, only linearizable functions are of interest and thus need be considered. If the intrinsic nonlinearity is signi®cant (which can be measured by eqn (9)), then f(b) cannot be nearly linearized and eqn (A.1) may not be nearly correct. As a result, con®dence region eqn (A.2) may not contain b with probability 1 ÿ a. Beale2 gave an approximate correction to obtain a conservative (1 ÿ a) ´ 100% con®dence region. With the correction, da2 is given by n N/ ps2 Fa p; n ÿ p; p 1 da2 1 nÿ1 n p 2 2 N/ ps2 Fa p; n ÿ p; p P 2 da 1 n ÿ pp A:6 where N/ is a theoretical measure of intrinsic nonlinearity obtained by letting sets bl be normally distributed, and then letting m ® 1 in (9)2 . (Actually, Beale2 made an approximation that Guttman and Meeter15 found to yield signi®cant error and Linssen20 corrected. The corrected version is used here.) In this study, we follow ^ / as given by Linssen's20 suggestion and replace N/ by M eqn (9). The signi®cance of the intrinsic nonlinearity is simply measured by the magnitude of the corrected da as compared to the uncorrected value. The same arguments as originally used to justify the use of eqn (A.5) to compute Schee-type intervals can be used to incorporate the correction eqn (A.6). Because the con®dence region is conservative, so are the Scheetype intervals. REFERENCES
Fig. 9. Relations between a con®dence interval (g(bL ), g(bU )) and possible true values g(b1 ) and g(b2 ) for g(b). S b ÿ S ^ b da2 is the con®dence region boundary.
1. Bates, D.M. & Watts, D.G., Relative curvature measures of nonlinearity. Journal of the Royal Statistical Society, Series B, 1980, 42(1), 1±25,. 2. Beale, E.M.L., Con®dence regions in non-linear estimation. Journal of the Royal Statistical Society, series B, 1960, 22(1), 41±76. 3. Beven, K. & Binley, A., The future of distributed models: Model calibration and uncertainty prediction. Hydrological Processes, 1992, 6, 279±298. 4. Brooks, R. J., Lerner, D. N. & Tobias, A. M., Determining the range of predictions of a groundwater model which arises from alternative calibrations. Water Resources Research, 1994, 30(11), 2993±3000. 5. Christensen, S., Hydrological model for the Tude aa catchment. Nordic Hydrology, 1994, 25(3), 145±166. 6. Christensen, S., Calibration data accuracy in a regionalscale groundwater model. Internal research report, Dept. of Earth Sciences, Aarhus University, Aarhus, 1995, 29 pp. 7. Christensen, S., On the strategy of estimating regionalscale transmissivity ®elds. Ground Water, 1997, 35(1), 131± 139.
Evaluation of con®dence intervals for a steady-state leaky aquifer model 8. Clarke, G.P.Y., Approximate con®dence limits for a parameter function in nonlinear regression. Journal of the American Statistical Association, 1987, 82(397), 221±230. 9. Cooley, R.L., Exact Schee-type con®dence intervals for output from groundwater ¯ow models, 1, Use of hydrogeologic information. Water Resources Research, 1993, 29(1), 17±33. 10. Cooley, R.L., Exact Schee-type con®dence intervals for output from groundwater ¯ow models, 2, Combined use of hydrogeologic information and calibration data. Water Resources Research, 1993, 29(1), 35±50. 11. Cooley, R.L. & Hill, M.C., A comparison of three Newton-like nonlinear least-squares methods for estimating parameters of ground-water ¯ow models. Computational Methods in Water Resources IX, Vol. 1, (ed. by T.F. Russell, R.E. Ewing, C.A. Brebbia, W.G. Gray & G.F. Pinder), 1992, 379±386. 12. Cooley, R.L. & Na, R.L., Regression modelling of ground-water ¯ow. U.S. Geological Survey Techniques of Water-Resources Investigations, 1990, book 3, ch. B4, 232 pp. 13. Cooley, R.L. & Vecchia, A.V., Calculation of nonlinear con®dence and prediction intervals for ground-water ¯ow models. Water Resources Bulletin, 1987, 23(4), 581±599. 14. Draper, N. & Smith, H., Applied regression analysis. 2nd ed., John Wiley, New York, 1981, 709 pp.
817
15. Guttman, I. & Meeter, D.A., On Beale's measures of nonlinearity. Technometrics, 1965, 7(4), 623±637. 16. Hamilton, D., & Wiens, D., Correction factors for F ratios in nonlinear regression. Biometrica, 1987, 74(2), 423±425. 17. Hill, M.C., Analysis of accuracy of approximate, simultaneous, nonlinear con®dence intervals on hydraulic heads in analytical and numerical test cases. Water Resources Research, 1989, 25(21), 177±190. 18. Hill, M.C., A computer program (MODFLOWP) for estimating parameters of a transient, three-dimensional ground-water ¯ow model using nonlinear regression. U.S. Geological Survey, Open-File Report 91-484, 1992, 358 pp. 19. Journel, A.G. & Huijbregts, C., Mining geostatistics. Academic Press, New York, 1978, 600 pp. 20. Linssen, H.N., Nonlinearity measures: A case study. Statistica Neerlandica, 1975, 29, 93±99. 21. Mood, A.M., Graybill, F.A. & Boes, D.C., Introduction to the theory of statistics. 3rd ed., McGraw-Hill, New York, 1974, 564 pp. 22. Seber, G.A.F. & Wild, C.J., Nonlinear regression. John Wiley, New York, 1989, 768 pp. 23. Vecchia, A.V. & Cooley, R.L., Simultaneous con®dence and prediction intervals for nonlinear regression models with application to a groundwater ¯ow model. Water Resources Research, 1987, 23(7), 1237±1250.