Use of PEST in Calibrating a Complex Groundwater Model

The Model

Figure 1 shows the domain of a finite difference groundwater model (MODFLOW) constructed to simulate a coastal aquifer system; a basemap is omitted from this figure to protect client confidentiality. All cells are 500m square; fixed head cells are highlighted.

Within the modeled area are a number of boreholes in which groundwater levels have been measured for the last 15 years. Locations of these bores are shown in Figure 1; an example hydrograph is shown in Figure 2. The pronounced seasonal variation of water table elevation is obvious from this figure.

Figure 1. MODFLOW grid, fixed heads cells and borehole
locations over study area.

Figure 2. Typical borehole hydrograph.

The groundwater model was built for two reasons:

• to predict water table variations resulting from possible future climatic scenarios, and

• to "process" groundwater data in order to draw inferences concerning recharge and evapotranspiration mechanisms operative within the catchment.

A simple recharge model (named LUMPREM for "LUMPed Parameter Recharge Model") was developed and applied to each of four zones defined within the model domain, zonation having been designed according to the clearing and subsequent re-vegetation history over different parts of the catchment. Figure 3 shows a schematic of the recharge model. (Note that modeled recharge and evapotranspiration rates from the soil moisture store bear a functional relationship to climate and current moisture content, the parameters of these relationships being partially determined through the model calibration process.) Figure 4 shows the four zones within the model domain to each of which a differently parameterized recharge model was attached, the differences in recharge resulting from differences in the amount and age of vegetation within each zone.

Figure 3. Schematic of the lumped parameter recharge model.
 

Figure 4. Zones within groundwater model domain to which
different recharge models were attached.

The model was comprised of a single, unconfined, MODFLOW "type 1" layer. However MODFLOW was altered in such a way that aquifer properties were allowed to vary with depth. Through this mechanism a single MODFLOW layer was used to simulate a regional two-layer system, it being obvious from borehole hydrographs that hydraulic conductivity is smaller in the lower part of the aquifer than it is in the upper part of the aquifer. (Anyone who has used the MODFLOW BCF or LPF packages will know of the instability and nonuniqueness problems to which it is prone if two MODFLOW layers are used to simulate a shallow aquifer system in which the water table regularly crosses the layer interface. By modifying MODFLOW in the manner described, the groundwater system in the study area could be simulated without the necessity of invoking two MODFLOW layers.)

Hydraulic property zonation is shown in Figure 5. In each of the 19 zones shown in this figure, estimates were required for upper and lower hydraulic conductivity as well as of upper and lower specific yield.

Figure 5. Zonation of upper and lower hydraulic
conductivity and specific yield.

Upon studying all available groundwater data, there were grounds for suspicion that vegetation met part of its evapotranspirational requirements by extracting water directly from the groundwater regime (or the capillary fringe directly above it). This phenomenon was simulated using the MODFLOW EVT package, requiring that evapotranspiration rates and extinction depths be estimated in each of the 16 zones shown in Figure 6, these zones having been selected according to vegetation type and age.

Figure 6. Zonation of parameters pertaining
to evapotranspiration.

Calibration Strategy

Great care must be taken in calibrating a model to resist attempts to estimate values for too many parameters. Where parameters are many, correlation between parameters is high and estimated values are unreliable. In the present case, even though borehole data was plentiful, some methodology was required to "regularize" (as the mathematicians say) the parameter estimation problem. As the discussion below demonstrates, two methods were employed. (Note that the work required to calibrate this model was carried out before the incorporation of advanced regularization functionality into PEST through its "regularization mode" option. Though the following discussion of regularization is sound, it could have been done even better using PEST, possibly with the help of pilot points for spatial parameter definition.)

The model was calibrated over a five year period from late 1991 to early 1997. In all, 120 parameter values required estimation. A total of 3542 water level reading were available from the boreholes scattered throughout the model domain. Using Parallel PEST, parameter values were adjusted until the fits between model-generated and borehole water levels were reduced to a minimum in the weighted least squares sense. Parallel PEST distributed model runs across three Pentium 150 machines running Microsoft WINDOWS 95. As the time required to run each model was about 10 minutes, use of Parallel PEST reduced optimization time from over 3 days to just one day; further reductions could have been achieved through using more machines.

The first regularization strategy involved the introduction of one extra "observation", added with appropriate weighting to the set of borehole water levels. PEST was informed that the "observed" value of recharge over the five-year simulation period within the large yellow zone of Figure 4 was 35% of rainfall over that time. "Tying down" recharge in this manner removed some of the correlation between model parameters. Without it, recharge model parameters together with hydraulic conductivity and specific yield values could all vary in harmony with little effect on model outcomes.

The second regularization strategy involved the insistence that groundwater model properties vary as little as possible over the model domain, compatible with allowing a faithful emulation of observed hydrographs. Hence a program was written to evaluate the standard deviation of upper and lower hydraulic conductivity and specific yield values over the model domain. An "observed value" of zero was supplied to PEST for each of these standard deviations and an appropriate weight was assigned. Hence PEST was asked to optimize parameter values in terms of reducing the discrepancies between simulated and observed borehole hydrographs under the condition that spatial parameter variation was kept to a minimum.

The "model" thus consisted of a four recharge models, a program to convert the results of each such recharge model to a form usable by the MODFLOW recharge package, MODFLOW itself, the spatial parameter standard deviation program, and the program used to interpolate MODFLOW heads to borehole locations (which is supplied with the PEST MODFLOW/MT3D Utilities); see Figure 7. A batch file was written to run each of these programs in the correct order; PEST ran this batch file whenever it needed to run the model.

Figure 7. The model as run by Parallel PEST.

Calibration Results

Figure 8 shows a typical fit obtained between measured and modeled water levels over the calibration period.

Figure 8. Modeled and measured borehole hydrographs
over calibration period.

The model was also run over a "verification period" consisting of the 6 years preceding the calibration period. Figure 9 shows a typical comparison of a measured and observed borehole hydrograph.

Figure 9. Modeled and measured borehole hydrographs
over verification period.

In the present circumstances, use of the groundwater model was not simply limited to its role as a predictor of groundwater response to different future climatic possibilities, for inferences could be drawn on the effect of vegetation on groundwater through PEST's parameterization of both the lumped parameter recharge models and the MODFLOW EVT package. In particular, it was possible to infer that:

• Clearing of vegetation results in a significant increase in recharge. However the greater the rainfall during a particular wet season, the smaller is the percentage increase. (This can be understood in terms of water moving quickly through a charged profile, or through soil macropores as upper subsurface conditions become saturated. Under these conditions downwardly-percolating water moves too fast for plant evapotranspiration to "grab" it, so that recharge becomes somewhat insensitive to the presence or absence of vegetation.)

• Recharge rates return to normal levels within about five years of revegetation.

• Native vegetation meets some of its evapotranspirational needs through direct extraction from the groundwater regime. However once an area has been revegetated, it takes at least 10 to 15 years for the new vegetation to have any impact on the groundwater regime through this mechanism.

• Within the study area, the upper 8 to 10 meters of the subsurface is on average 10 times more hydraulically conductive than the underlying subsurface.

Conclusions

Without a sophisticated model-independent, nonlinear parameter estimator such as PEST, calibration of the complex recharge/groundwater model incorporating two important regularizing schemes would have been impossible. The spatial complexity of the model parameterization allows accurate simulation of groundwater behavior over all parts of its regime; yet the fact that parameter values have been constrained to vary as little as necessary to achieve this simulation accuracy removes the possibility of large spatial parameter "oscillation" to which unregularized parameter estimation is prone.

Nonlinear parameter estimation techniques are indispensable for model calibrating and data interpretation. PEST has brought with it a quantum leap in the applicability of this methodology to the day-to-day problems that beset modelers from all field of science and engineering by facilitating deployment of this methodology in conjunction with existing models, and by allowing the development of innovative and powerful parameterization strategies such as have been demonstrated here.