Systematic Parameter Estimation Strategy for Refining the Birkenes Model

S Y S T E M A T I C P A R A M E T E R E S T I M A T I O N S T R A T E G Y F O R R E F I N I N G T H E B I R K E N E S M O D E L

J . W . Delleur

School of Civil Engineering, Purdue University, West Lafayette, IN 47907, USA

F i - J o h n C h a n g

Department of Agricultural Engineering, National Taiwan University, Taipei, Taiwan, ROC

A B S T R A C T A systematic parameter estimation strategy

which includes a regionalized sensitivity analysis and an a u t o m a t i c p a r a m e t e r calibration technique was developed to improve p a r a m e t e r estimation. Results obtained using the Birkenes model of watershed acidification show t h a t the p a r a m e t e r s of the model can then be systematically estimated, leading to better model performance.

I N T R O D U C T I O N

The Birkenes model was developed, using a simple two-reservoir hydrogeochemical conceptual model, for a forested catchment in southernmost Norway for predicting the acidification of soils and freshwaters (Christophersen et al., 1982). Although the model s t r u c t u r e is r a t h e r simple and general, its application to a certain watershed requires t h a t the values of its p a r a m e t e r s be made specific for the given situation. Since the accuracy of the model performance depends upon the adequacy of the estimation of its p a r a m e t e r s , the availability of a reliable parameter estimation technique will dominate its usefulness in future applications.

Two types of estimation procedures are commonly used: the traditional, m a n u a l procedures and the newer, automatic procedures. The m a n u a l estimation is a trial and error procedure; it strongly depends on the user's experience with the model and the watershed which is being modeled. This method is obviously tedious and time consuming. In contrast, the automatic procedure has been found to be more efficient when an a p p r o p r i a t e systematic optimization technique can be developed.

However, automatic procedures have some inherent difficulties to overcome: for example, the interdependence of some of the parameters and the uncertainty in finding the global optimum in the presence of possible local optima for the objective function surface. In an a t t e m p t to overcome these difficulties, a systematic p a r a m e t e r estimation strategy is proposed which consists of a Regionalized Sensitivity Analysis (RSA) (Hornberger and Spear, 1981) followed by an Automatic P a r a m e t e r Calibration Technique (APCT). The RSA eliminates the insensitive parameters and establishes the a p p r o p r i a t e range of values for each remaining parameter. Consequently, the optimal solution is more likely to be obtained by an A P C T method. This joint s t r a t e g y is implemented for the Birkenes model, to simplify its use and to improve its performance.

D A T A U S E D

The observed d a t a sets (météorologie, hydrological, and streamflow chemical d a t a ) for the three years, 1977 through 1979, in a subcatchment known as "Inflow No. 4" of H a r p Lake, 8 km northeast of Huntsville, Ontario,

164 J.W. Delleur and F.J. Chang

C a n a d a , were obtained from the Ontario Ministry of Environment, Dorset Research Centre, C a n a d a . The first year d a t a were used for RSA and A P C T calibration, and the last two years of d a t a were used for verification in this study. The subcatchment has been described by Seip et al. (1985) and in papers referenced therein.

R E G I O N A L I Z E D S E N S I T I V I T Y A N A L Y S I S ( R S A )

The procedure of this method is given by the following steps: Step 1: Select the parameters to be tested.

Step 2: Define the design ranges of the selected p a r a m e t e r s .

Step 3: For each selected parameter, generate a series of independent r a n d o m numbers with a uniform distribution within the design range. T h e series typically numbered 300.

Step 4: Randomly select a sample from each series to form a p a r a m e t e r set and run the model using this particular set of p a r a m e t e r values. Whether a p a r a m e t e r set is acceptable or unacceptable is based on a subjective R value, which represents the model improvement over the mean of the observed values. If the R value obtained from the simulation is greater than this subjective R value, the result is acceptable, otherwise, the result is unacceptable. The simulation R value is calculated as follows:

E(x-xi)

- Eft-xO

R = i = i — - i=i (i)

E(x-x0

i = l

where n is the number of observations; x; are observed values, i—1 to n; x is their mean; and x, are the predicted values.

Step 5: For each parameter, compare the distribution of the p a r a m e t e r values associated with the acceptable results to the distribution for the unacceptable results. The Kolmogorov-Smirnov (K-S) test with a specified level of significance (e.g. 5%) is used to classify whether or not the distributions for the acceptable and unacceptable results are significantly different. If the two functions are not statistically different, then the parameter is insensitive for the simulation of t h a t particular set of results; otherwise the p a r a m e t e r is sensitive. For the insensitive p a r a m e t e r s , specific values are assigned within the previous design ranges. For the sensitive p a r a m e t e r s , the range of values is reduced to a favorable range which, in general, corresponds to the steepest portion of the acceptable cumulative distribution.

A U T O M A T I C P A R A M E T E R C A L I B R A T I O N T E C H N I Q U E ( A P C T )

The procedure is summarized in the following steps: Step 1: Select the parameters to be optimized. Step 2: Define a range of values for each parameter. Step 3: Give an initial set of the p a r a m e t e r values. Step 4: Define the objective function (Maximize R).

Step 5: Optimize each p a r a m e t e r successively. In this study, the Golden Section method is used to find the p a r a m e t e r value giving the maximum R value. After all parameters have been optimized, one optimization cycle is complete. The value of the objective function, (R), is stored at the end of

the cycle.

Step 6: Repeat the optimization cycle (Step 5) until the R value does not increase by more t h a n a small amount, for example, 0.001.

Flow charts of the RSA and A P C T procedures can be found in C h a n g (1988).

A P P L I C A T I O N O F R S A M E T H O D T O H Y D R O L O G I C A L P O R T I O N O F B I R K E N E S M O D E L

Since the hydrological flow p a t h s are apparently critical to the u n d e r s t a n d i n g of the changes in the chemical characteristics of surface and subsurface waters (Driscoll and Newton, 1985), the hydrological portion of the Birkenes model ( Figure 1) is fitted first.

AMIN I-Ea P

t . I

If

F i g u r e 1 Schematic of t h e Birkenes hydrologie submodel.

Seven p a r a m e t e r s are selected; they are: A, B, Ak, Bk, Amin, Bmax, and Bmin. A and B are t h e upper and lower reservoir volumes (mm), respectively. Ak and Bk are discharge r a t e constants (1/day). Amin, Bmin, and Bmax are reservoir threshold p a r a m e t e r s . The design ranges of the p a r a m e t e r values are given in Table 1. The original values of the

166 J.W. Delleur and F.J. Chang

T a b l e 1. Results of K-S test for the Birkenes model (Hydrologie portion)

p a r a m e t e r A B Ak Bk Amin Bmax Bmin original value 13.0 112 0.7 0.04 15.0 160 100 design range 7.5-26 56-224 0.4-1.0 0.02-0.08 7.5-30 80-320 50-200 P(1) 0.790 0.000 0.0005 0.063 0.958 0.000 0.206 sensi-tive NO YES YES YES NO YES NO favorable range 13.0 70-140 0.48-0.72 0.025-0.06 15.0 225-325 100.0 (1): The p value is defined as the level of significance

at which the K-S test would just pass. parameters were obtained from Seip et a l , (1986).

Here, a subjective value of R = 0 . 4 is used, and yields approximately the same numbers of acceptable and unacceptable responses. The cumulative distributions of the acceptable and unacceptable parameters are shown in Figure 2. The K-S test a t the significance level 5 % was used to classify whether or not the two distributions are significantly different. T h e results are shown in Table 1. The two distributions for Ak, B, and Bmax are significantly different; the two distributions for A, Amin, and Bmin are not significantly different. The p value of Bk is slightly greater t h a n 5%, so it is also classified as a sensitive parameter. For each of the insensitive parameters, a specific value within previously defined ranges was assigned. For the sensitive parameters Ak, B, Bk and Bmax, the ranges were reduced to the favorable ranges obtained from Figure 2 and are listed in Table 1.

A P P L I C A T I O N OF A P C T T O T H E H Y D R O L O G I C A L P O R T I O N OF T H E B I R K E N E S M O D E L

Four cases were considered. As shown in Table 2, cases 1 and 2 included seven p a r a m e t e r s . The p a r a m e t e r ranges of case 1 were the same as in the Monte-Carlo simulation in the previous section, while case 2 was set to correspond to the favorable range of the RSA. T h e optimal R value of case 2 was the same as t h a t of case 1 but was obtained in less t h a n half the number of iterations (model simulations).

Cases 3 and 4 were set to include only the four sensitive p a r a m e t e r s with their favorable ranges as identified by the RSA in the previous section. Both cases, each with different initial values, yielded the same final optimal values of the p a r a m e t e r s . T h u s , the optimal solution was likely a global one ( or very close to the global solution) for this particular p a r a m e t e r set. Moreover, even though their optimal R values are slightly less t h a n in case

1, fewer iterations were needed to reach the optimal solution for cases 3 and 4. The reuslts d e m o n s t r a t e t h a t , because the number of iterations is greatly reduced, the joint RSA-AFCT strategy is of great worth. This is a critical consideration as each iteration requires a complete simulation of the integrated acidification model.

Parameter Estimation for Birkenes Model 16 7 l 0 0 o \o o \ ^ V \ \ N . \ ^% •> Tj " O » o o « 5 O cq (X) d faO o -a in a .2 -fa » 3 £ 1 'C -fa » .2 3 •3 a> -fa »

a,

168 J.W. Delleur and F.J. Chang

T a b l e 2 Results of A P C T for the hydrological portion

case 1 2 3 4 source initial range final l e n g t h '1' optimal value initial range final length optimal value initial range final length optimal value initia! range final length optimal value A 13 7.5-26 1.0 26.0 7.5 7.5-15 1.0 15.0 13.0 _(2) (2) 13.0 13.0 13.0 B 112 56-224 4.0 89.87 70 70-140 4.0 119.57 70 0 70-140 4.0 96.738 140.0 70-140 4.0 96.738 Ak 0.7 0.4-1.0 0.005 0.4 0.4 0.4-0.7 0.005 0.4 0.48 0.48-0.72 0.005 0.48 0.72 0.48-0.72 0.005 0.48 Bk 0.04 0.02-0.08 0.004 0.0463 0.025 0.25-0.06 0.004 0.0466 0.026 0.026-0.06 0.004 0.0371 0.06 0.026-0.06 0.004 0.0371 Amin 15.0 7.5-30 1.0 10.78 7.5 7.5-15 1.0 10.365 15.0 15.0 15.0 15.0 Bmax 160.0 80-320 4.0 320 225 225-320 4.0 320 225.0 225-325 4.0 325.0 325.0 225-325 4.0 325.0 Bmin 100.0 50-200 4.0 137.54 90.0 90-115 4.0 155 100.0 100.0 100.0 100.0 N( 3 ) R 682 0.6433 310 0.6453 73 0.6113 68 0.6113

(l); Minimum uncertainty interval for golden section search. (2): Not calibrated. (3): Number of iterations

T a b l e 3 Results of 3 years simulation using the optimal p a r a m e t e r values

case 0 1 2 o u A 13.0 26.0 15.0 13.0 B 112 89.873 119.6 96.738 Ak 0.7 0.4 0.40 0.48 Bk 0.04 0.046 .047 .037 Amin 15.0 10.783 10.365 15.0 Bmax 160 320 320 325 Bmin 100 137.6 155 100 R 0.4943 0.677 0.6783 0.6753

case 0, used by Seip et al., (1986), gave R = 0 . 4 9 , while the p a r a m e t e r values of cases 1 through 3, obtained by the A P C T method, gave R values greater than 0.67. Comparing the R values in cases 1 through 3 with the values in case 0, a significant improvement (about 18%) in the model performance was shown. This demonstrates t h a t the automatic p a r a m e t e r calibration technique not only provides an efficient way of calibrating the p a r a m e t e r s , but also significantly improves the model performance in the streamflow prediction.

S Y S T E M A T I C P A R A M E T E R E S T I M A T I O N S T R A T E G Y F O R T H E C H E M I C A L P O R T I O N OF B I R K E N E S M O D E L

In this section, the model performance of the streamflow sulfate concentration was tested. The regionalized sensitivity analysis was used to identify the sulfate parameters of the Birkenes model. Six parameters which correlate with streamflow sulfate concentration were chosen. These were: FA, FB, EQCA, EQCB, CKB, and CBLIK. FA and F B are the total amounts of water-soluble sulfate in the upper and lower soil reservoirs, respectively. E Q C A is a proportionality constant. CKB is a r a t e constant. CBLIK is an assumed equilibrium concentration in the lower reservoir. E Q C B is a proportionality constant.

In the calculation of the criterion R, the variables (x, x, x;) in eq. 1 were used for the streamflow sulfate concentrations. Three subjective R values were used to classify the results as acceptable or unacceptable. The results are shown in Table 4.

T a b l e 4 The K-S test results of the Birkenes model (Chemical portion)

parameter FA FB EQCA EQCB CKB CBLIK original value 235000 1025000 1200.0 5000.0 0.02 200.0 design range 1762500-293750 768750-1281250 900-1500 3750-6250 0.015-0.025 150-250 accepted responses unaccepted responses R = 0 . 3 P 0.589 0.058 0.165 0.486 0.746 0.000 159 141 R=0.35 P 0.287 0.049 0.205 0.314 0.282 0.000 144 156 R = 0 . 4 P 0.21 0.12 0.267 0.38 0.74 0.000 sensi-tive NO NO NO NO NO YES 124 176

T h e three R values produced similar results of p a r a m e t e r classifications. This indicates t h a t the criterion is robust. The parameter CBLIK shows very small p values in all cases, so it was classified as a sensitive parameter. The other p a r a m e t e r s were classified as insensitive parameters since their p values were larger t h a n 5%.

The A P C T was then used to calibrate the sensitive parameter CBLIK. T h e A P C T procedure yielded an optimal R value of 0.591. The model was simulated for the next two years for verification as shown in Figure 3. The streamflow sulfate concentrations from model predictions fitted the observed values reasonably well ( R = 0 . 4 3 ) . Details of the results can be found in C h a n g (1988).

C O N C L U S I O N S

The a u t o m a t i c p a r a m e t e r calibration technique can be implemented for the Birkenes model. The A P C T procedure becomes the main program and the model is modified to a subroutine. The p a r a m e t e r s of the refined Birkenes model can then be calibrated automatically without m a n u a l adjustment and

17 0 J.W. Delleur and F.J. Chang ~i/oa n

a better model performance obtained. In an a t t e m p t to overcome the problems of inherent nonlinear properties in the model s t r u c t u r e and multiple local optima of the objective function surface, it is recommended t h a t the regionalized sensitivity analysis be executed to eliminate the insensitive p a r a m e t e r s and establish the favorable simulation ranges of the sensitive p a r a m e t e r s before the A P C T is performed.

A C K N O W L E D G E M E N T The authors wish to t h a n k P.J. Dillon and Ed de Grosbois, Dorset Research Centre, Ontario Ministry of the Environment, C a n a d a , for making the Harp Lake d a t a available.

R E F E R E N C E S

Chang, Fi-John, (1988) Refinement of Hydrogeochemical Models of the Ecological Impact of Acid Deposition. Ph.D Thesis, P u r d u e University. Christophersen, N., Seip, H.M. and Wright, R.F., (1982) A model for Stream

W a t e r Chemistry a t Birkenes, Norway. Wat. Resour. Res., 18(4), 977-996.

Driscoll, C.T. and Newton, R.M., (1985) Chemical Characteristics Adirondack Lakes. Environ. Sci. Technol., 19, 1018-1024.

Hornberger, G.M. and Spear, R.C., (1981) An Approach of the Preliminary Analysis of Environmental Systems. /. of Environ. Management, 12,

7-18.

Seip, H.M., Seip, R., Dillon, P.J. and de Grosbois, Ed, (1985) Model of Sulfate Concentration in a Small Stream in the Harp Lake Catchment, Ontario, Can. J. Fish Aquat. Sci., 42, 927-937.