- Email: [email protected]
Application of four watershed acidification models to Batchawana Watershed, Canada
Environmental Pollution 77 (1992) 243-252
Application of four watershed acidification models to Batchawana Watershed, Canada W. G. Booty, A. G. Bobba, D. C. L. Lam & D. S.
Jeff'ties
Environment Canada, National Water Research Institute, Rivers Research Branch, 867 Lakeshore Road, Burlington, Ontario, Canada L7R 4A6
Four watershed acidification models (TMWAM, ETD, ILWAS, and RAINS) are reviewed and a comparison of model performance is presented for a common watershed. The models have been used to simulate the dynamics of water quantity and quality at Batchawana Watershed, Canada, a sub-basin of the Turkey Lakes Watershed. The computed results are compared with observed data for a fouryear period (Jan. 1981-Dec. 1984). The models exhibit a significant range in the ability to simulate the daily, monthly and seasonal changes present in the observed data. Monthly watershed outflows and lake chemistry predictions are compared to observed data. pH and ANC are the only two chemical parameters common to all four models. Coefficient of efficiency (E), linear (r) and rank (R) correlation coefficients, and regression slope (s) are used to compare the goodness of fit of the simulated with the observed data. The ILWAS, TMWAM and RAINS models performed very well in predicting the monthly flows, with values of r and R of approximately 0.98. The ETD model also showed strong correlations with linear (r) and rank (R) correlation coefficients of 0.896 and 0.892, respectively. The results of the analyses showed that TMWAM provided the best simulation of pH (E = 0-264, r = 0-648), which is slightly better than ETD (E = 0.240, r = 0-549), and much better than ILWAS (E = -2.965, r = 0.293), and RAINS (E---4-004, r = 0-473). ETD was found to be superior in predicting ANC (E = 0.608, r = 0-781) as compared to TMWAM (E -- 0.340, r = 0.598), ILWAS (E = 0.275, r = 0-442), and RAINS (E = - 1.048, r -- 0.356). The TMWAM model adequately simulated SO4 over the four-year period (E = 0.423, r -- 0.682) but the ETD ( E = - 0 . 9 0 4 , r = 0.274), ILWAS (E=-4-314, r = 0-488), and RAINS (E = -6-479, r = 0.126) models all performed poorer than the benchmark model (mean observed value).
INTRODUCTION
Enhanced Trickle-Down (ETD) model (Schnoor et al., 1984; Nikolaidas et ai., 1988), the M A G I C model (Cosby et al., 1985a) and the Regional Acidification Information and Simulation Model System (RAINS) (Kamari, 1985). The rate of reaching the equilibrium or steady-state concentration for a given acid load differs from lake to lake and if the loading is reduced the recovery times will also vary. Theoretically, the watershed response times range from decades to centuries (Arp & Ramnarine, 1983; Cosby et al., 1985b; Davis & Goldstein, 1988; Ormerod et al., 1988). The I L W A S and M A G I C models are being used in the N A P A P Integrated Assessment in the US and the R A I N S model is being used for a similar study in Europe. Validation of these models on a wide range o f lake-watershed systems is an important process which is used to evaluate the expected accuracy and uncertainty of the predictions generated by the models. As part of the 1990 Canadian
A watershed modeling approach is required to be able to predict the effects of varying levels of acidic precipitation on aquatic and terrestrial systems. A number of process-oriented mathematical models have been developed that simulate the temporal hydrological and chemical responses o f surface waters. Some examples of these mechanistic models include the Integrated Lake Watershed Acidification Study (ILWAS) model (Chen et al., 1983), T u r k e y - M e r s e y Watershed Acidification Model ( T M W A M ) (Bobba & Lam, 1988; L a m et al., 1988), the Birkenes model (Christophersen et al., 1982; Beck et al., 1987), the Watershed Acidification Model (WAM) (Booty, 1983; Booty & Kramer, 1984), the Environ. Pollut. 0269-7491/92/$05.00 © 1992 Elsevier Science Publishers Ltd, England. Printed in Great Britain 243
244
W. G. Booty, A. G. Bobba, D. C. L. Lam, D. S. Jeffries
Long-Range Transport of Air Pollutants and Acid D e p o s i t i o n A s s e s s m e n t ( R M C C , 1990), it was determ i n e d t h a t it w o u l d be useful to e x a m i n e the a c c u r a c y o f a n u m b e r o f these m o d e l s by c o m p a r i n g t h e m to a c o m m o n d a t a base. A s a case study, the following m o d e l s ( T M W A M , E T D , I L W A S a n d R A I N S ) were a p p l i e d to the T u r k e y L a k e s W a t e r s h e d (Fig. 1) in C a n a d a for s h o r t - t e r m (four years) analysis a n d the M A G I C m o d e l for l o n g - t e r m (50-100 years) analysis. This p a p e r presents the results o f the s h o r t - t e r m analysis.
simplifications o f reality, a n d uncertainties inherent in the structure o f the m o d e l s as well as uncertainties in the characteristics o f the w a t e r s h e d s being m o d e l e d prevent the m o d e l s f r o m being used to m a k e a b s o l u t e p r e d i c t i o n s o f w a t e r s h e d responses to a t m o s p h e r i c inputs. H o w e v e r , b y c o m p a r i n g the responses o f different m o d e l s a p p l i e d to a c o m m o n watershed, a feeling for the influence o f the m o d e l design on p r e d i c t i n g w a t e r s h e d responses can be obtained.
Conceptual aspects MODEL COMPARISONS T a b l e 1 presents a c o m p a r i s o n o f the m a j o r processes a n d the levels o f r e s o l u t i o n in the f o u r models, for the versions t h a t were c u r r e n t in 1989. T h e m o d e l s were originally d e v e l o p e d to satisfy different objectives for different g e o g r a p h i c regions. C o n s e q u e n t l y , the m o d e l s reflect these differences in the c o m p l e x i t y o f process f o r m u l a t i o n s used, a n d levels o f spatial a n d t e m p o r a l a g g r e g a t i o n . All o f the models, b y definition, are
Table 1. Major processes included in the four models (Modified from Jenne et ai., 1989) Watershed/lake processes
Models ILWAS ETD T M W A M RAINS
Canopy Dry scavenging Exudation/neutralization
+ +
-
-
-
+ + + + + + + + + +
+ + + + + + -
+ + + + + + + -
+ + + + + + + +
+ + + + + + + +
+ + + + -
+ + + + + + -
+ + + + + -
+ + + + + + +
+ + + + -
+ -
+ -
Hydrology Interception Evapotranspiration Vertical infiltration Root zone flow Lateral interflow Groundwater Snowmelt Streamflow Lake stratification Lake routing
Terrestrial biogeochernical Ion exchange Mineral weathering Carbonate reactions Sulfate reactions Organic acid reactions Nitrogen reactions Aluminum reactions Decomposition/respiration
Lake~stream biogeochemical Ion exchange Mineral weathering Carbonate reactions Sulfur reactions Nitrogen reactions Organic acid reactions Aluminum reactions
There are three basic w a y s o f designing a w a t e r s h e d acidification model. T h e first m e t h o d uses the simplest m o d e l t h a t is c a p a b l e o f representing the variables o f interest. T h e second m e t h o d is to d e v e l o p a m o d e l o f i n t e r m e d i a t e c o m p l e x i t y that is n o r m a l l y designed for a limited range o f watersheds. T h e third m e t h o d is to design a general m o d e l t h a t is c o m p r e h e n s i v e e n o u g h to be a p p l i e d to a n y watershed. T h e T M W A M , E T D , a n d R A I N S m o d e l s were originally d e v e l o p e d using the second m e t h o d while the I L W A S m o d e l was d e v e l o p e d using the third m e t h o d . The T M W A M consists o f four m a j o r s u b m o d e l s that deal with h y d r o l o g y , sulfate, c a t i o n s a n d D O C . It is a l u m p e d p a r a m e t e r model, as are the o t h e r three m o d e l s considered in this paper. The u n d e r l y i n g h y p o t h e s i s is t h a t the w a t e r s h e d h y d r o l o g y is largely responsible for the d y n a m i c v a r i a t i o n s a n d t h a t the w a t e r a n d soil chemistry follow simple equilibrium conditions. T M W A M consists o f a h y d r o l o g i c a l m o d e l ( B o b b a & L a m , 1988) which was interfaced to a m o d i f i e d version o f the geochemical s u b m o d e l o f the Birkenes m o d e l ( C h r i s t o p h e r s e n et al., 1982). The details o f the T M W A M m o d e l have been described earlier in the literature ( B o b b a & L a m , 1988, 1989, 1990; L a m e t al., 1988, 1989). The E T D m o d e l is a c o m p a r t m e n t a l i z e d mass b a l a n c e a p p r o a c h to estimating the c o n c e n t r a t i o n s o f chloride, sulfate, a l k a l i n i t y a n d p H in sensitive w a t e r s h e d s receiving acid deposition. It consists o f a m a s s balance for a l k a l i n i t y (acid neutralizing c a p a c i t y ) in several c o m p a r t m e n t s ( u p p e r soil, u n s a t u r a t e d zone, s a t u r a t e d zone, b e d r o c k aquifer, a n d lake c o m p a r t m e n t s ) . The R A I N S m o d e l consists o f four modules. T h e m o d e l c o m p u t e s the lake acidity d i s t r i b u t i o n s for the lakes within a given lake region as a function o f w a t e r s h e d characteristics a n d local deposition. D e t a i l s o f the m o d e l are presented by K a m a r i (1985). E a c h w a t e r s h e d is divided c o n c e p t u a l l y into four spatial sectors: s n o w p a c k , u p p e r soil layer, lower soil layer, a n d lake volume. Different m o d u l e s o f the m o d e l c o m p u t e the h y d r o l o g y a n d the flux o f ions c o n t r i b u t i n g to the acidity a n d a l k a l i n i t y o f the lake water. The I L W A S m o d e l was d e v e l o p e d as a process o r i e n t e d m o d e l t h a t c o u l d be used as a tool to e v a l u a t e the interconnections a n d interactions o f a t m o s p h e r i c t e r r e s t r i a l - a q u a t i c systems as p a r t o f the overall process
Application of watershed acidification models of surface water acidification. The model was originally developed for application to watersheds in the Adirondack region but because it is based on general principles of continuity and conservation of mass, it is generally applicable to most lake-watershed systems. The model simulates the biogeochemical processes that occur along the pathways which precipitation follows as it passes through the forest canopy, soil horizons, streams, bogs, and lakes. ILWAS is basically a lumped parameter model that attempts to minimize the heterogeneity of watershed parameters by dividing the watershed into horizontally 'homogeneous' sub-basins, stream and lake segments. Each sub-basin is segmented vertically into compartments representing canopy, snowpack and soil layers. The lake is also segmented vertically for the simulation of heat and chemical fluxes.
Atmospheric processes Each of the models requires data on precipitation quantity and air temperature. TMWAM, ILWAS, and ETD models require daily meteorological data whereas the RAINS model uses monthly data. Precipitation can be broken down into rain or snow using air temperatures, which are also used to determine the snowmelt rates in all the models. Only the ILWAS model uses daily air temperatures to calculate the lakewater and soil temperatures.
Hydrologic processes The four models utilize different levels of aggregation of hydrologic processes and different levels of complexity in the functions used to simulate the processes. The ILWAS model employs the greatest level of detail and considers all of the hydrologic processes listed in Table 1. The other three models are intermediate in complexity (with respect to state-of-the-art hydrologic modeling). Only the ILWAS model considers canopy interception. In the other three models it is assumed that canopy storage is not a significant component of the hydrologic budget or that its affect on precipitation volume reaching the soil can be empirically represented. Evapotranspiration is determined as a function of the potential evapotranspiration in all four models. Each model uses a different empirical function, but most are based on air temperature, relative humidity, latitude, and wind speed. All of the functions that are used contain terms that are based on monthly averaged data. In the TMWAM model, the soils through which water flows are divided into three compartments, called the upper, lower and groundwater reservoirs. In the ETD model, the subsurface hydrologic system is broken into the upper soil, the unsaturated zone, the saturated zone, and the bedrock aquifer compartments. The RAINS model considers two soil layers (upper and lower). The hydrologic module routes precipitation into surface
245
runoff, base flow and interlayer flows. Since the flowpaths determine which materials water contacts in its passage through the watershed and the length of time the water is in contact with these materials, predicting the chemistry of the water in these systems should be highly dependent on the accurate prediction of the flow routing. The accuracy with which predictions of hydrologic flowpaths can be made is limited by the spatial variability of the hydraulic characteristics throughout the watershed. Intuitively, the most accurate prediction should be generated by the most highly parameterized process representations. However, the more densely parameterized codes tend to contain variables which are not easily measured or which are based on limited aggregated observations. In the Batchawana Watershed, approximately 35% of the precipitation occurs in the form of snow, and spring snowmelt is a major component of the annual water budget. Snowpack volume is simulated by all four models. Each model uses an empirical expression that relates the snowmelt rate to air temperature. The ILWAS, TMWAM, and ETD models also consider the decrease of the snowpack through sublimation. Only the ILWAS model considers the hydraulics of flow within a stream channel. However, due to the fact that the streams within Batchawana Watershed are first order with poorly defined channel geometries, the stream hydrology option was not used. Batchawana Lake (north and south basins) exhibits dimictic stratification, with partial to complete mixing in the spring and fully mixed conditions in the autumn. Only the ILWAS and RAINS models consider stratification within a lake. The ILWAS code allows for up to 80 horizontal layers while the RAINS code allows for two. The ILWAS model also simulates the temperature profile within a lake and the formation of an ice layer. As snowmelt normally occurs before the ice cover on a lake thaws and breaks up, snowmelt is commonly observed flowing over or just under the ice cover on the lake. This short-circuiting of meltwater at Batchawana Lake has been observed (Jeffries & Semkin, 1983) and has been shown to significantly affect lakewater chemistry.
Biogeochemical processes Terrestrial All of the models incorporate functions to represent ion exchange, mineral weathering, carbonate, and sulfate reactions. The ILWAS and TMWAM models also consider organic acid reactions. In the Batchawana Watershed, the concentrations of organic acids are negligible and consequently the lack of these species being considered in the other models was not found to be a serious weakness for this particular application. Only the ILWAS model attempts to simulate the nitrogen cycle. Studies of ion balances in the watershed (Jeffries et al., 1988) indicated that NH4 and NO3 contribute 10%
246
W. G. Booty, A. G. Bobba, D. C. L. Lam, D. S. Jeffries
to the total sum of cations and anions, respectively. It was also found for the Batchawana Watershed that the molar ratio of anions to cations was closer to unity if only the major ions were considered rather than all of the ions. Difficulties in measuring forest uptake rates of the N compounds was a key reason for excluding them from the TMWAM model during its development and application to the Turkey Lakes Watershed ( L a m et al., 1988). Aluminum species are considered in all of the models except ETD. In Batchawana Watershed, surface and subsurface water pH values normally range between 5.5 and 7-0. In this pH range, the contribution of aluminum species to total ANC is very small to negligible. Decomposition and respiration processes are only emulated in the ILWAS model. Root respiration and litter decay generate CO2(g) in the soil zone, which is either lost to the atmosphere or dissolved in the soil solution. Advective transport of CO2(g) is calculated from the percolation rates of the fluid phase in the soil layers and from changes in moisture content. The other models use fixed CO2(g) partial pressures for the different soil horizons simulated.
Aquatic Ion exchange, mineral weathering, and sulfur reactions are not considered within stream or lakewaters in the TMWAM or RAINs models, and only within the lake for the ETD model. Due to the short lake-residence times for the two Batchawana Lake basins (1.3 and 0-3 years), in-lake processes, such as sulfate reduction, are not considered to be significant processes affecting lake ion budgets (Jeffries et al., 1988). In all of the model codes, the carbonate reactions are considered equilibrium processes. The ILWAS and ETD models use similar functions to describe the effect of temperature on the equilibrium constants. The TMWAM and RAINS models do not correct the equilibrium constants for changing soil and lakewater temperatures.
Input data requirements The ILWAS model is by far the most complex of the four models and it is also the most input data intensive. Considerable effort was required to obtain the necessary physical, chemical, and biological data for Batchawana Watershed that are required for the numerous input files. Indeed, more time was spent generating the input files for the ILWAS model than for the other three models combined. This is due to the fact that the ILWAS model allows a much greater degree of detail. For soils in Batchawana Watershed, it was determined that the soils could be divided into five layers which are chemically unique. The other models only allow two or three subsurface layers to be defined. The ILWAS model allows the user to use up to 10 subsurface layers and provides for up to 20 subcatchments and associated
stream segments. The other model codes do not provide an option for the model user to break the watershed into subcatchments. In the authors' application of the ILWAS model to the Batchawana Watershed, the watershed was divided into two subcatchments, one for Batchawana Lake N. and one for Batchawana Lake S. (Fig. 1).
Numerical methods In the ILWAS model, the time step is variable with a maximum time step of one day. Very short time steps are not necessary for model stability or for convergence. However, many of the input data, such as precipitation chemistry, are supplied in monthly averaged form and many of the model process coefficients are also supplied as monthly values. The ETD and TMWAM models also use a one-day time step to calculate the hydrologic flow routing and rate-limited chemical processes. The ETD model uses a predictor-corrector integration technique to solve the first-order differential equations. A Runge-Kutta algorithm is used to initiate calculations at the beginning of each time step. This method is less efficient than the standard Eulerian numerical integration techniques used in ILWAS and results in the model requiring a proportionately greater length of computing time than the other models, considering the size of the source code. The TMWAM model uses analytical techniques to solve the first-order differential flow equations. In the geochemical submodel a Newton iterative method is used to solve a fourth-order equation with respect to H ÷.
Calibration procedure Hydrologic calibration For all the models, the hydrologic components are calibrated before attempting to calibrate the chemistry. A total of nine parameters were adjusted during the calibration of the ILWAS model hydrologic module for the Batchawana Watershed. For the ETD model, 13 parameters were optimized during the calibration of the hydrologic portion of the model. For TMWAM and RAINS, six parameters were optimized in order to carry out calibrations of the hydrologic modules. Chemical calibration Calibration of the chemical modules of the models was carried out for all of the master chemical variables that are used in each model. However, emphasis was placed on pH, alkalinity, and sulfate, as these were the three common parameters. In calibrating the chemical module of the ILWAS model, 143 parameters were optimized. However, many of these were determined using a separate program supplied by the model developers, called SOILCALC. In comparison, for the ETD model, 17 parameters needed to be optimized and both the
Application of watershed acidification models 64026 ,
I I \
47 °
84025 . I
247
64024 . I
/
/
/
f
84 ° 23' I
f~'="--~X
\
Batchawana
Lake N.~
.47°04 '
• \
-...(
I /
Batchat
S1
\
Lake S. /
/
/ Lake
i=S4
art
f \
\
\
, \
I
I /I
.47003 e
/
/
Approximate Watershed Boundary
)" 250
/ f °"t°'i° U O"e"c I
0
250 500 750 1000
Metres
I-">~,,,.,,=' -4,,.,
.47002 '
• Stream gauging and sampling station
/~SA
~
J ,4128'
~,~25'
84~24 '
84~23 '
Fig. 1. Location map of the Turkey Lakes Watershed showing catchment boundary and stream gauging/sampling stations. TMWAM and RAINS models have only five parameters that had to be optimized.
APPLICATION OF MODELS TO THE TURKEY LAKES (BATCHAWANA LAKE SUB-BASIN) The ILWAS, RAINS, ETD, and TMWAM models were applied to the Batchawana sub-basin in order to evaluate and compare the accuracy of the models in predicting several common hydrological and chemical parameters. Batchawana sub-basin is the headwater subbasin of the Turkey Lakes Watershed and is approximately 2.06 km 2 in surface area and receives approximately 1212 mm of precipitation annually. A detailed description of the Turkey Lakes Watershed and its subbasins has been given in Jeffries and Semkin (1982). The ILWAS, ETD and TMWAM models operate on a daily time step but the RAINS model operates on a monthly time step. Also, as not all of the models are capable of simulating all of the major cations and
anions, pH, alkalinity, and sulfate, which are common to all the models, are discussed for comparison purposes• All of the models were applied for the period 1 Jan. 1981-31 Dec. 1984 for data measured at station SI, shown in Fig. 1. The models were initially calibrated for the period 1 Jan. 1981-31 Dec. 1982. They were then run for the final three years of data without further calibration• The models were run using the same basic physical and chemical input constants and parameter values as well as meteorological input data.. Every effort was made to calibrate the different models to the best degree of accuracy, and, where available, the model developer's procedure for model calibration was followed. The models were all calibrated by trial and error rather than by an automatic numerical optimization procedure. The model parameters were adjusted to minimize the errors between the observed and model calculated outputs, based on the principle of least squares. Additional statistics are used to evaluate model performance. Linear and rank correlation coefficients and regression slopes based on a statistical comparison
248
W. G. Booty, A. G. Bobba, D. C. L. Lain, D. S. Jeffries
of computed and observed values are also used to indicate model performance. In the rank correlation, the correlation coefficient (r) and the regression slope (s) are derived from the ordered ranks of both the measured and computed data and are good indicators of episodic coincidences. A correlation coefficient and slope close to 1.0 are indicative of good model performance (Bobba & Lam, 1988). The coefficient of efficiency (E) (Bobba et al., 1988) is a measure of the improvement on the variance not explainable by the mean value, but explainable by the particular model which produces the computed values (qc). The coefficient quantifies the relative improvement of the predictions over a benchmark model (i.e. the mean observed value). The coefficient of efficiency (E) (Nash & Sutcliffe, 1970; Bobba & Lam, 1988) is defined as:
E~
predicting monthly outflows from the watershed. The hydrologic submodel used in the ETD model is intermediate in complexity between the ILWAS and TMWAM and RAINS models. However, for this particular application, it shows the largest deviation from the measured values for both the monthly and cumulative flows, as reflected by all four model performance statistics. Considering each model performance statistic separately, it can be seen that, with respect to the linear correlation coefficient, the RAINS and TMWAM models performed the best with values of r -- 0-987 and 0.986, respectively, followed closely by the ILWAS model with r -- 0.967, and the ETD model was the least precise with r = 0.896. The rank correlation coefficient indicates how well the models simulate the timing and magnitude of the peaks of flows as compared to the linear correlation coefficient. In this case, the TMWAM, ILWAS, and RAINS models perform equally well, with R = 0-988, 0.985, and 0.984, respectively. Considering the coefficient of efficiency, the RAINS and TMWAM models perform the best with respect to the mean-asmodel, with E = 0.970 and 0.961, respectively, followed closely by the ILWAS model with R = 0.915, and finally the ETD model with R -- 0.762. The models show a significant range in their abilities to simulate the monthly and seasonal changes represented by the observed data for outlet chemistry. The model predicted versus observed outlet stream pH, alkalinity, and sulfate data are presented in Figs 3-5, respectively, and the associated statistical data on model performance with respect to these chemical parameters are shown in Table 2. Visually, it can be seen in Figs 3 and 4 that the RAINS model tends to overpredict the fluctuations of observed pH and alkalinity. The ETD model tends to underpredict the fluctuations of observed
[~(qo,--~lo)2--~(qo,-2]-q¢,) i=1
i=1
E (qo,-
(1) 2
i=1
where qo is the mean of the observed values. Any positive value of E can be interpreted as improvement over the mean-as-model, the closer to + 1 the better. The model simulated versus observed monthly watershed outflow data are presented in Fig. 2 and the associated model performance satatistical data are presented in Table 2. The TMWAM and RAINS models show very similar simulated flows and have almost identical statistical performance values. The ILWAS model, which has the most sophisticated hydrologic submodel, performed equally as well as the TMWAM and RAINS models, except with respect to its coefficient of efficiency, in 600
500
400 ~1<
"
300
,T
200 ,-.-.
,oo 0
///
),.. ........I, -
-
1981 •
Obs.Data
1982 TMV~M
1983 ----
ETD
: ......
1984 RAINS
ILWAS
Fig. 2. Comparison of model simulated versus observed monthly outflow data for Batchawana Lake Watershed at station S1 for the period 1 Jan. 1981-31 Dec. 1984.
Application of watershed acidification models
249
Table 2. Statistical evaluation of model performance for (a) the calibration period (1 Jan. 1981-31 Dec. 1982); (b) period 1 Jan. 1983-31 Dec. 1984; and (c) period 1 Jan. 1981-31 Dec. 1984 Model
Linear correlation coefficient (r)
Rank correlation coefficient (R)
ILWAS TMWAM ETD RAINS
0.967 0.986 0-896 0.987
0.985 0.988 0.892 0.984
ILWAS TMWAM ETD RAINS
0-293 0.438 0.549 0.473
0.613 0-713 0-477 0-177
ILWAS TMWAM ETD RAINS
0-442 0.598 0.781 0.356
0.511 0.533 0.687 0.407
ILWAS TMWAM ETD RAINS
0.488 0.682 0.274 0.126
0.610 0.680 0-233 0.153
Slope (s)
Coefficient of efficiency (E)
Outflow 1.100 0.912 1-181 1.065
0.915 0.961 0.762 0.970
-6-772 0.648 1.616 0.142
-2.965 0.264 0-240 -4.004
0.441 0.851 0.959 0.281
0-275 0.340 0.608 -1.048
0.334 0.819 0.444 -0.172
-4-314 0.423 -0.904 -6-479
pH
A/k
SO,
Values of r and R between 0.8 and 1.0 are considered a strong correlation. A value of s between 0.9 and 1-1 is considered to be a strong correlation. The coefficient of efficiency (E) is a measure of the improvement on the variance. Any positive value of E can be interpreted as improvement over the meanas-model, the closer to + 1 the better. pH but does the best job of simulating the observed alkalinities. The T M W A M model appears to do the best job of simulating the observed pH data and is second to the E T D model in simulating the measured alkalinity data. As seen in Table 2, the T M W A M model sulfate predictions have the highest correlation with the ob-
served data (r = 0-682, R = 0.680, s = 0.819). It is also the only model that has a positive coefficient of efficiency ( E = 0.423). In this case, the other three models' predictions of sulfate concentration in the outlet stream cannot be considered as improvement over the mean-asmodel. In Fig. 5 the observed versus model-predicted data are presented for the concentrations of sulfate at
6.5
6
-r" 5.5
5 ¸
4.5
1981
Oba.Data
1982
- -
TMWAM
1983
---
ETD
....... R A I N S
1984
- -
ILWAS
Fig. 3. Comparison of model simulated versus observed monthly averaged outflow pH values for Batchawana Lake Watershed at station SI for the period 1 Jan. 1981-31 Dec. 1984.
W. G. Booty, A. G. Bobba, D. C. L. Lain, D. S. Jeffries
250 150
,j
/i
100
g~
co
//i
50
<
0
I 1981
Obe.Data
Fig. 4.
I
1982
TMWAM
1983
-
--
ETD
1984
RAINS
- -
ILWAS
Comparison of model simulated versus observed monthly averaged outflow alkalinities for Batchawana Lake Watershed at station S1 for the period 1 Jan. 1981-31 Dec. 1984.
station S I. TMWAM predictions most closely simulate both the magnitude and the timing of the observed sulfate stream concentration peaks, as confirmed by the rank correlation coefficient of 0-680 (Table 2). The ILWAS model also performs well in predicting the timing of the peaks in stream sulfate concentration, with R--0.610, as shown in Table 2. However, the magnitudes of the peaks, especially in the spring of 1983 and 1984, are significantly different from the observed values. As was the case for pH, the ETD model predictions for SO4 show the smallest fluctuations of magnitude over the period of study. However, the timing of the peaks do not coincide very well with those of the observed data, as seen by the low rank correlation coefficient (R = 0-233).
The RAINS model underpredicts the stream sulfate concentrations over the entire period of study. It also fails to predict the timing of the peaks, as verified by the very low rank correlation coefficient shown in Table 2 (R = 0.153).
DISCUSSION The models use different definitions for ANC, as listed in Table 3. This may result in the values being different from those obtained using the strict Gran titration definition for the proton reference level. However, for the pH range (5.5-6.7) observed in the lake outflows, the
200
150
,j
2 100 (o c~ n -I 09
5O
I
1981 Oba.Dat8
Fig. 5.
I
1982 - -
TMWAM
1983 ---
ETD
........ R A I N S
1984 - -
ILWAS
Comparison of model simulated versus observed monthly averaged outflow sulfate concentrations for Batchawana Lake Watershed at station S1 for the period 1 Jan. 1981-31 Dec. 1984.
Application of watershed acidification models Table 3. Model definitions of acid neutralizing capacity (ANC) (Modified from Jenne et al., 1989)
Model
ANC definition
ILWAS
ANC -- [HCO~] + 2[CO 2-] + [OH-] + [H2R'"-] + 2[HR '''2-] + 3[R''a-] + JR'-] + [AIOH 2+] + 2[AI(OH)~] + 3[AI(OH)3°] + 4[AI(OH)g] + 3[AIR''~] + lAIR '2+] + 2[AI(R')~] + 3[AI(R')~ - [H +]
ETD
ANC = [HCO3] + 2[CO~ ] + [OH-] - [H +] + [R' ]
251
structure are responsible for the different levels of performance shown by the models. Many models are calibrated and verified over similar short periods of observation and then are used to make future predictions for longer periods of time (i.e. 10, 20, 50 years). It is clearly advantageous to evaluate the relative uncertainties in model predictions using multiple model predictions for measured data before making long-term forecasts of model responses to changes in such inputs as acid loadings.
TMWAM ANC = ([H ÷] + [Na +] + 2[Ca2+] + 2[Mg2+]) - ([HCO;] + 2[SO;])
RAINS
ANC = [HCO;] - [H +] - [AP÷]
Square brackets indicate molar units of concentration, and R', R", and R'" represent monoprotic, diprotic, and triprotic organic acids, respectively.
differences in the ANC values calculated by ILWAS and T M W A M should be comparable to the analytically defined A N C values. The ANC values presented for the measured data were determined using a Metrohm model 636 Titroprocessor (Jeffries et al., 1988) which detects the titration inflection point from the second derivative of the curve. Rigorous checking of this methodology has shown that the results obtained are equivalent to those determined using the FI Gran function. No correction factors have been applied to any of the model-predicted ANC values used in the comparisons presented in this paper. As it has been shown in numerous studies (Johnson et al., 1969; Jeffries & Semkin, 1982; Bobba et al., 1986; Rustad et al., 1986), watershed stream chemistry is strongly correlated with stream hydrology. It might be assumed that, if the hydrological submodel of the ETD model was improved so that its performance was more similar to those shown in Table 2 for the other three models, its ability to predict stream chemistry parameters such as pH, alkalinity, and sulfate may also be improved. The results of these comparisons of model performance for Batchawana Lake Watershed outflow are presented only as an example of how well the models perform for this one site, at one point, and based on monthly averaged data. If we were to apply the models to another period of observations, and/or to different watershed parameters, the relative performance of the models may be expected to be different than those presented in this study. However, the main point is to show that basic watershed parameters, such as basin outflow, outflow pH, alkalinity, and sulfate are simulated to different degrees of accuracy by the four different models considered in this study. A great deal more effort is still required in order to fully test and validate the models, and more importantly, to help to understand which key processes or differences in model
ACKNOWLEDGMENTS The assistance of R. Semkin in the acquisition of watershed data is gratefully acknowledged. The authors are also very grateful to H. M. Seip, N. Christophersen, S. A. Gherini, R. K. Munson, M. M. Lang, J. L. Schnoor, N. P. Nickolaidas and J. Kamari for providing the various model codes and for useful discussions.
REFERENCES
Arp, P. A. & Ramnarine, S. (1983). Verifying the effects of acid precipitation on soil leachates--A comparison between published records and model predictions. Ecol. Modelling, 19, 119-38. Beck, M. B., Drummond, D., Kleissen, F. M., Langan, S. J., Wheater, H. S. & Whitehead, P. G. (1987). Surface water acidification: The Birkenes data revisited. In Systems Analysis in Water Quality Management, ed. M. B. Beck, IAWPRC/Pergamon, New York, pp. 133-50. Bobba, A. G. & Lam, D. C. L. (1988). Application of a hydrological model to the Turkey Lakes watershed. Can. J. Fish. Aquat. Sci., 45, 81-7. Bobba, A. G. & Lam, D. C. L. (1989). Application of a hydrological model to acidified watersheds: A study on Mersey River and Moosepit Brook, Nova Scotia. Water, Air & Soil Poll., 46, 261-75. Bobba, A. G. & Lam, D. C. L. (1990). Hydrological modelling of acidified Canadian watersheds. Ecol. Modelling, 50, 532. Bobba, A. G., Lam, D. C. L. & De Grosbois, E. (1986). Application of hydrological model to acidified watersheds: A study on Harp Lake catchments. Hydrol. Sci. Technol. Short Pap., 2, 37-44. Booty, W. G. (1983). Watershed acidification model and the soil acid neutralization capacity concept, PhD dissertation, McMaster University, Hamilton, Ontario. Booty, W. G. & Kramer, J. R. (1984). Sensitivity analysis of a watershed acidification model. Philos. Trans. R. Soc. London, B305, 441-9. Chen, C. W., Gherini, S. A., Hodson, R. J. M. & Dean, J. D. (1983). The ILWAS study: Model principles and applications procedures. EPRI EA-3221, Project 1109-5, Vol. 1, Electric Power Research Institute, Palo Alto, CA. Christophersen, N., Seip, H. M. & Wright, R. F., (1982). A model for streamwater chemistry at Birkenes, Norway. Water Resour. Res., 18, 977-96. Cosby, B. J., Hornberger, G. M., Galloway, J. N. & Wright, R. F. (1985a). Modelling the effects of acid deposition: Assessment of a lumped parameter model of soil water and stream chemistry. Water Resour. Res., 21, 51-63.
252
W. G. Booty, A. G. Bobba, D. C. L. Lam, D. S. Jeffries
Cosby, B. J., Hornberger, G. M., Galloway, J. N. & Wright, R. F. (1985b). Time scales of catchment acidification: A quantitative model for estimating freshwater acidification. Environ. Sct Technol., 19, 1144-9. Davis, J. E. & Goldstein, R. A. (1988). Simulated response of an acidic Adirondack lake-watershed to various liming mitigation strategies. Water Resour. Res., 24, 525-32. Jeffries, D. S. & Semkin, R. G. (1982). Basin description and information pertinent to mass balance studies of the Turkey Lakes watershed. Turkey Lakes Watershed Rep. TLW-8201, Can. Centre Inland Waters, Burlington, Ontario. Jeffries, D. S. & Semkin, R. G. (1983). Changes in snowpack, stream, and lake chemistry during snowmelt in the Turkey Lakes Watershed. VDI-Berichte, 500, 377-86. Jeffries, D. S., Semkin, R. G., Neureuther, R. & Seymour, M. (1988). Ion mass budgets for lakes in the Turkey Lakes watershed, June 1981-May, 1983. Can. J. Fish. Aquat. Sci., 45, 47-58. Jenne, E. A. et al. (1989). An evaluation and analysis of three dynamic watershed acidification codes (MAGIC, ETD, and ILWAS). US Environmental Protection Agency, Washington, DC, PNL-6687, UC-11. Johnson, N. M., Likens, G. E., Bormann, F. H., Fisher, D. W. & Pierce, R. S. (1969). A working model for the variation of stream water chemistry at the Hubbard Brook Experimental Forest, New Hampshire. Water Resour. Res., 5, 1353-63. Kamari, J. (1985). A model for analyzing lake water acidification on a regional scale. Report no. CP-85-48. International Institute for Applied Systems Analysis, 2361 Laxenberg, Austria.
Lam, D. C. L., Bobba, A. G., Jeffries, D. S. & Craig, D. (1988). Modelling stream chemistry for the Turkey Lakes Watershed: Comparison with 1981-84 data. Can. J. Fish. Aquat. Sci., 45, 72-80. Lam, D. C. L., Bobba, A. G., Bourbonniere, R. A., Howell, G. D. & Thompson, M. E. (1989). Modelling organic and inorganic acidity in two Nova Scotia rivers. Water, Air & Soil Poll., 46, 277-87. Nash, J. E. & Sutcliffe, J. V. (1970). Riverflow forecasting through conceptual models. 1. A discussion of principles. J. Hydrol., 10, 282-90. Nikolaidas, N. P., Rajaram, H., Schnoor, J. L. & Geogakakos, K. P. (1988). A generalized soft water acidification model. Water Resour. Res., 24, 1983-96. Ormerod, S. J., Weatherly, N. S. & Whitehead, P. G. (1988). Temporal patterns of ecological change during the acidification and recovery of streams throughout Wales according to a hydrochemical model. In Regional Acidification Models, ed. J. Kamarai. Springer-Verlag, New York, pp. 167-84. RMCC (Research Monitoring Coordinating Committee) (1990). Canadian Long-Range Transport of Air Pollutants and Acid Deposition Assessment Report. Volumes !-8, Environment Canada, Ottawa. Rustad, S., Christophersen, N., Seip, H. M. & Dillon, P. J. (1986). Model for streamwater chemistry of a tributary to Harp Lake, Ontario. Can. J. Fish. Aquat. Sci., 43, 625-33. Schnoor, J. L., Palmer, W. D., Jr & Glass, G. E. (1984). Modeling impacts of acid precipitation for northeastern Minnesota. In Modeling of Total Acid Precipitation Impact, ed. J. L. Schnoor. Butterworth, Boston, MA, pp. 155-73.
Application of four watershed acidification models to Batchawana Watershed, Canada W. G. Booty, A. G. Bobba, D. C. L. Lam & D. S.
Jeff'ties
Environment Canada, National Water Research Institute, Rivers Research Branch, 867 Lakeshore Road, Burlington, Ontario, Canada L7R 4A6
Four watershed acidification models (TMWAM, ETD, ILWAS, and RAINS) are reviewed and a comparison of model performance is presented for a common watershed. The models have been used to simulate the dynamics of water quantity and quality at Batchawana Watershed, Canada, a sub-basin of the Turkey Lakes Watershed. The computed results are compared with observed data for a fouryear period (Jan. 1981-Dec. 1984). The models exhibit a significant range in the ability to simulate the daily, monthly and seasonal changes present in the observed data. Monthly watershed outflows and lake chemistry predictions are compared to observed data. pH and ANC are the only two chemical parameters common to all four models. Coefficient of efficiency (E), linear (r) and rank (R) correlation coefficients, and regression slope (s) are used to compare the goodness of fit of the simulated with the observed data. The ILWAS, TMWAM and RAINS models performed very well in predicting the monthly flows, with values of r and R of approximately 0.98. The ETD model also showed strong correlations with linear (r) and rank (R) correlation coefficients of 0.896 and 0.892, respectively. The results of the analyses showed that TMWAM provided the best simulation of pH (E = 0-264, r = 0-648), which is slightly better than ETD (E = 0.240, r = 0-549), and much better than ILWAS (E = -2.965, r = 0.293), and RAINS (E---4-004, r = 0-473). ETD was found to be superior in predicting ANC (E = 0.608, r = 0-781) as compared to TMWAM (E -- 0.340, r = 0.598), ILWAS (E = 0.275, r = 0-442), and RAINS (E = - 1.048, r -- 0.356). The TMWAM model adequately simulated SO4 over the four-year period (E = 0.423, r -- 0.682) but the ETD ( E = - 0 . 9 0 4 , r = 0.274), ILWAS (E=-4-314, r = 0-488), and RAINS (E = -6-479, r = 0.126) models all performed poorer than the benchmark model (mean observed value).
INTRODUCTION
Enhanced Trickle-Down (ETD) model (Schnoor et al., 1984; Nikolaidas et ai., 1988), the M A G I C model (Cosby et al., 1985a) and the Regional Acidification Information and Simulation Model System (RAINS) (Kamari, 1985). The rate of reaching the equilibrium or steady-state concentration for a given acid load differs from lake to lake and if the loading is reduced the recovery times will also vary. Theoretically, the watershed response times range from decades to centuries (Arp & Ramnarine, 1983; Cosby et al., 1985b; Davis & Goldstein, 1988; Ormerod et al., 1988). The I L W A S and M A G I C models are being used in the N A P A P Integrated Assessment in the US and the R A I N S model is being used for a similar study in Europe. Validation of these models on a wide range o f lake-watershed systems is an important process which is used to evaluate the expected accuracy and uncertainty of the predictions generated by the models. As part of the 1990 Canadian
A watershed modeling approach is required to be able to predict the effects of varying levels of acidic precipitation on aquatic and terrestrial systems. A number of process-oriented mathematical models have been developed that simulate the temporal hydrological and chemical responses o f surface waters. Some examples of these mechanistic models include the Integrated Lake Watershed Acidification Study (ILWAS) model (Chen et al., 1983), T u r k e y - M e r s e y Watershed Acidification Model ( T M W A M ) (Bobba & Lam, 1988; L a m et al., 1988), the Birkenes model (Christophersen et al., 1982; Beck et al., 1987), the Watershed Acidification Model (WAM) (Booty, 1983; Booty & Kramer, 1984), the Environ. Pollut. 0269-7491/92/$05.00 © 1992 Elsevier Science Publishers Ltd, England. Printed in Great Britain 243
244
W. G. Booty, A. G. Bobba, D. C. L. Lam, D. S. Jeffries
Long-Range Transport of Air Pollutants and Acid D e p o s i t i o n A s s e s s m e n t ( R M C C , 1990), it was determ i n e d t h a t it w o u l d be useful to e x a m i n e the a c c u r a c y o f a n u m b e r o f these m o d e l s by c o m p a r i n g t h e m to a c o m m o n d a t a base. A s a case study, the following m o d e l s ( T M W A M , E T D , I L W A S a n d R A I N S ) were a p p l i e d to the T u r k e y L a k e s W a t e r s h e d (Fig. 1) in C a n a d a for s h o r t - t e r m (four years) analysis a n d the M A G I C m o d e l for l o n g - t e r m (50-100 years) analysis. This p a p e r presents the results o f the s h o r t - t e r m analysis.
simplifications o f reality, a n d uncertainties inherent in the structure o f the m o d e l s as well as uncertainties in the characteristics o f the w a t e r s h e d s being m o d e l e d prevent the m o d e l s f r o m being used to m a k e a b s o l u t e p r e d i c t i o n s o f w a t e r s h e d responses to a t m o s p h e r i c inputs. H o w e v e r , b y c o m p a r i n g the responses o f different m o d e l s a p p l i e d to a c o m m o n watershed, a feeling for the influence o f the m o d e l design on p r e d i c t i n g w a t e r s h e d responses can be obtained.
Conceptual aspects MODEL COMPARISONS T a b l e 1 presents a c o m p a r i s o n o f the m a j o r processes a n d the levels o f r e s o l u t i o n in the f o u r models, for the versions t h a t were c u r r e n t in 1989. T h e m o d e l s were originally d e v e l o p e d to satisfy different objectives for different g e o g r a p h i c regions. C o n s e q u e n t l y , the m o d e l s reflect these differences in the c o m p l e x i t y o f process f o r m u l a t i o n s used, a n d levels o f spatial a n d t e m p o r a l a g g r e g a t i o n . All o f the models, b y definition, are
Table 1. Major processes included in the four models (Modified from Jenne et ai., 1989) Watershed/lake processes
Models ILWAS ETD T M W A M RAINS
Canopy Dry scavenging Exudation/neutralization
+ +
-
-
-
+ + + + + + + + + +
+ + + + + + -
+ + + + + + + -
+ + + + + + + +
+ + + + + + + +
+ + + + -
+ + + + + + -
+ + + + + -
+ + + + + + +
+ + + + -
+ -
+ -
Hydrology Interception Evapotranspiration Vertical infiltration Root zone flow Lateral interflow Groundwater Snowmelt Streamflow Lake stratification Lake routing
Terrestrial biogeochernical Ion exchange Mineral weathering Carbonate reactions Sulfate reactions Organic acid reactions Nitrogen reactions Aluminum reactions Decomposition/respiration
Lake~stream biogeochemical Ion exchange Mineral weathering Carbonate reactions Sulfur reactions Nitrogen reactions Organic acid reactions Aluminum reactions
There are three basic w a y s o f designing a w a t e r s h e d acidification model. T h e first m e t h o d uses the simplest m o d e l t h a t is c a p a b l e o f representing the variables o f interest. T h e second m e t h o d is to d e v e l o p a m o d e l o f i n t e r m e d i a t e c o m p l e x i t y that is n o r m a l l y designed for a limited range o f watersheds. T h e third m e t h o d is to design a general m o d e l t h a t is c o m p r e h e n s i v e e n o u g h to be a p p l i e d to a n y watershed. T h e T M W A M , E T D , a n d R A I N S m o d e l s were originally d e v e l o p e d using the second m e t h o d while the I L W A S m o d e l was d e v e l o p e d using the third m e t h o d . The T M W A M consists o f four m a j o r s u b m o d e l s that deal with h y d r o l o g y , sulfate, c a t i o n s a n d D O C . It is a l u m p e d p a r a m e t e r model, as are the o t h e r three m o d e l s considered in this paper. The u n d e r l y i n g h y p o t h e s i s is t h a t the w a t e r s h e d h y d r o l o g y is largely responsible for the d y n a m i c v a r i a t i o n s a n d t h a t the w a t e r a n d soil chemistry follow simple equilibrium conditions. T M W A M consists o f a h y d r o l o g i c a l m o d e l ( B o b b a & L a m , 1988) which was interfaced to a m o d i f i e d version o f the geochemical s u b m o d e l o f the Birkenes m o d e l ( C h r i s t o p h e r s e n et al., 1982). The details o f the T M W A M m o d e l have been described earlier in the literature ( B o b b a & L a m , 1988, 1989, 1990; L a m e t al., 1988, 1989). The E T D m o d e l is a c o m p a r t m e n t a l i z e d mass b a l a n c e a p p r o a c h to estimating the c o n c e n t r a t i o n s o f chloride, sulfate, a l k a l i n i t y a n d p H in sensitive w a t e r s h e d s receiving acid deposition. It consists o f a m a s s balance for a l k a l i n i t y (acid neutralizing c a p a c i t y ) in several c o m p a r t m e n t s ( u p p e r soil, u n s a t u r a t e d zone, s a t u r a t e d zone, b e d r o c k aquifer, a n d lake c o m p a r t m e n t s ) . The R A I N S m o d e l consists o f four modules. T h e m o d e l c o m p u t e s the lake acidity d i s t r i b u t i o n s for the lakes within a given lake region as a function o f w a t e r s h e d characteristics a n d local deposition. D e t a i l s o f the m o d e l are presented by K a m a r i (1985). E a c h w a t e r s h e d is divided c o n c e p t u a l l y into four spatial sectors: s n o w p a c k , u p p e r soil layer, lower soil layer, a n d lake volume. Different m o d u l e s o f the m o d e l c o m p u t e the h y d r o l o g y a n d the flux o f ions c o n t r i b u t i n g to the acidity a n d a l k a l i n i t y o f the lake water. The I L W A S m o d e l was d e v e l o p e d as a process o r i e n t e d m o d e l t h a t c o u l d be used as a tool to e v a l u a t e the interconnections a n d interactions o f a t m o s p h e r i c t e r r e s t r i a l - a q u a t i c systems as p a r t o f the overall process
Application of watershed acidification models of surface water acidification. The model was originally developed for application to watersheds in the Adirondack region but because it is based on general principles of continuity and conservation of mass, it is generally applicable to most lake-watershed systems. The model simulates the biogeochemical processes that occur along the pathways which precipitation follows as it passes through the forest canopy, soil horizons, streams, bogs, and lakes. ILWAS is basically a lumped parameter model that attempts to minimize the heterogeneity of watershed parameters by dividing the watershed into horizontally 'homogeneous' sub-basins, stream and lake segments. Each sub-basin is segmented vertically into compartments representing canopy, snowpack and soil layers. The lake is also segmented vertically for the simulation of heat and chemical fluxes.
Atmospheric processes Each of the models requires data on precipitation quantity and air temperature. TMWAM, ILWAS, and ETD models require daily meteorological data whereas the RAINS model uses monthly data. Precipitation can be broken down into rain or snow using air temperatures, which are also used to determine the snowmelt rates in all the models. Only the ILWAS model uses daily air temperatures to calculate the lakewater and soil temperatures.
Hydrologic processes The four models utilize different levels of aggregation of hydrologic processes and different levels of complexity in the functions used to simulate the processes. The ILWAS model employs the greatest level of detail and considers all of the hydrologic processes listed in Table 1. The other three models are intermediate in complexity (with respect to state-of-the-art hydrologic modeling). Only the ILWAS model considers canopy interception. In the other three models it is assumed that canopy storage is not a significant component of the hydrologic budget or that its affect on precipitation volume reaching the soil can be empirically represented. Evapotranspiration is determined as a function of the potential evapotranspiration in all four models. Each model uses a different empirical function, but most are based on air temperature, relative humidity, latitude, and wind speed. All of the functions that are used contain terms that are based on monthly averaged data. In the TMWAM model, the soils through which water flows are divided into three compartments, called the upper, lower and groundwater reservoirs. In the ETD model, the subsurface hydrologic system is broken into the upper soil, the unsaturated zone, the saturated zone, and the bedrock aquifer compartments. The RAINS model considers two soil layers (upper and lower). The hydrologic module routes precipitation into surface
245
runoff, base flow and interlayer flows. Since the flowpaths determine which materials water contacts in its passage through the watershed and the length of time the water is in contact with these materials, predicting the chemistry of the water in these systems should be highly dependent on the accurate prediction of the flow routing. The accuracy with which predictions of hydrologic flowpaths can be made is limited by the spatial variability of the hydraulic characteristics throughout the watershed. Intuitively, the most accurate prediction should be generated by the most highly parameterized process representations. However, the more densely parameterized codes tend to contain variables which are not easily measured or which are based on limited aggregated observations. In the Batchawana Watershed, approximately 35% of the precipitation occurs in the form of snow, and spring snowmelt is a major component of the annual water budget. Snowpack volume is simulated by all four models. Each model uses an empirical expression that relates the snowmelt rate to air temperature. The ILWAS, TMWAM, and ETD models also consider the decrease of the snowpack through sublimation. Only the ILWAS model considers the hydraulics of flow within a stream channel. However, due to the fact that the streams within Batchawana Watershed are first order with poorly defined channel geometries, the stream hydrology option was not used. Batchawana Lake (north and south basins) exhibits dimictic stratification, with partial to complete mixing in the spring and fully mixed conditions in the autumn. Only the ILWAS and RAINS models consider stratification within a lake. The ILWAS code allows for up to 80 horizontal layers while the RAINS code allows for two. The ILWAS model also simulates the temperature profile within a lake and the formation of an ice layer. As snowmelt normally occurs before the ice cover on a lake thaws and breaks up, snowmelt is commonly observed flowing over or just under the ice cover on the lake. This short-circuiting of meltwater at Batchawana Lake has been observed (Jeffries & Semkin, 1983) and has been shown to significantly affect lakewater chemistry.
Biogeochemical processes Terrestrial All of the models incorporate functions to represent ion exchange, mineral weathering, carbonate, and sulfate reactions. The ILWAS and TMWAM models also consider organic acid reactions. In the Batchawana Watershed, the concentrations of organic acids are negligible and consequently the lack of these species being considered in the other models was not found to be a serious weakness for this particular application. Only the ILWAS model attempts to simulate the nitrogen cycle. Studies of ion balances in the watershed (Jeffries et al., 1988) indicated that NH4 and NO3 contribute 10%
246
W. G. Booty, A. G. Bobba, D. C. L. Lam, D. S. Jeffries
to the total sum of cations and anions, respectively. It was also found for the Batchawana Watershed that the molar ratio of anions to cations was closer to unity if only the major ions were considered rather than all of the ions. Difficulties in measuring forest uptake rates of the N compounds was a key reason for excluding them from the TMWAM model during its development and application to the Turkey Lakes Watershed ( L a m et al., 1988). Aluminum species are considered in all of the models except ETD. In Batchawana Watershed, surface and subsurface water pH values normally range between 5.5 and 7-0. In this pH range, the contribution of aluminum species to total ANC is very small to negligible. Decomposition and respiration processes are only emulated in the ILWAS model. Root respiration and litter decay generate CO2(g) in the soil zone, which is either lost to the atmosphere or dissolved in the soil solution. Advective transport of CO2(g) is calculated from the percolation rates of the fluid phase in the soil layers and from changes in moisture content. The other models use fixed CO2(g) partial pressures for the different soil horizons simulated.
Aquatic Ion exchange, mineral weathering, and sulfur reactions are not considered within stream or lakewaters in the TMWAM or RAINs models, and only within the lake for the ETD model. Due to the short lake-residence times for the two Batchawana Lake basins (1.3 and 0-3 years), in-lake processes, such as sulfate reduction, are not considered to be significant processes affecting lake ion budgets (Jeffries et al., 1988). In all of the model codes, the carbonate reactions are considered equilibrium processes. The ILWAS and ETD models use similar functions to describe the effect of temperature on the equilibrium constants. The TMWAM and RAINS models do not correct the equilibrium constants for changing soil and lakewater temperatures.
Input data requirements The ILWAS model is by far the most complex of the four models and it is also the most input data intensive. Considerable effort was required to obtain the necessary physical, chemical, and biological data for Batchawana Watershed that are required for the numerous input files. Indeed, more time was spent generating the input files for the ILWAS model than for the other three models combined. This is due to the fact that the ILWAS model allows a much greater degree of detail. For soils in Batchawana Watershed, it was determined that the soils could be divided into five layers which are chemically unique. The other models only allow two or three subsurface layers to be defined. The ILWAS model allows the user to use up to 10 subsurface layers and provides for up to 20 subcatchments and associated
stream segments. The other model codes do not provide an option for the model user to break the watershed into subcatchments. In the authors' application of the ILWAS model to the Batchawana Watershed, the watershed was divided into two subcatchments, one for Batchawana Lake N. and one for Batchawana Lake S. (Fig. 1).
Numerical methods In the ILWAS model, the time step is variable with a maximum time step of one day. Very short time steps are not necessary for model stability or for convergence. However, many of the input data, such as precipitation chemistry, are supplied in monthly averaged form and many of the model process coefficients are also supplied as monthly values. The ETD and TMWAM models also use a one-day time step to calculate the hydrologic flow routing and rate-limited chemical processes. The ETD model uses a predictor-corrector integration technique to solve the first-order differential equations. A Runge-Kutta algorithm is used to initiate calculations at the beginning of each time step. This method is less efficient than the standard Eulerian numerical integration techniques used in ILWAS and results in the model requiring a proportionately greater length of computing time than the other models, considering the size of the source code. The TMWAM model uses analytical techniques to solve the first-order differential flow equations. In the geochemical submodel a Newton iterative method is used to solve a fourth-order equation with respect to H ÷.
Calibration procedure Hydrologic calibration For all the models, the hydrologic components are calibrated before attempting to calibrate the chemistry. A total of nine parameters were adjusted during the calibration of the ILWAS model hydrologic module for the Batchawana Watershed. For the ETD model, 13 parameters were optimized during the calibration of the hydrologic portion of the model. For TMWAM and RAINS, six parameters were optimized in order to carry out calibrations of the hydrologic modules. Chemical calibration Calibration of the chemical modules of the models was carried out for all of the master chemical variables that are used in each model. However, emphasis was placed on pH, alkalinity, and sulfate, as these were the three common parameters. In calibrating the chemical module of the ILWAS model, 143 parameters were optimized. However, many of these were determined using a separate program supplied by the model developers, called SOILCALC. In comparison, for the ETD model, 17 parameters needed to be optimized and both the
Application of watershed acidification models 64026 ,
I I \
47 °
84025 . I
247
64024 . I
/
/
/
f
84 ° 23' I
f~'="--~X
\
Batchawana
Lake N.~
.47°04 '
• \
-...(
I /
Batchat
S1
\
Lake S. /
/
/ Lake
i=S4
art
f \
\
\
, \
I
I /I
.47003 e
/
/
Approximate Watershed Boundary
)" 250
/ f °"t°'i° U O"e"c I
0
250 500 750 1000
Metres
I-">~,,,.,,=' -4,,.,
.47002 '
• Stream gauging and sampling station
/~SA
~
J ,4128'
~,~25'
84~24 '
84~23 '
Fig. 1. Location map of the Turkey Lakes Watershed showing catchment boundary and stream gauging/sampling stations. TMWAM and RAINS models have only five parameters that had to be optimized.
APPLICATION OF MODELS TO THE TURKEY LAKES (BATCHAWANA LAKE SUB-BASIN) The ILWAS, RAINS, ETD, and TMWAM models were applied to the Batchawana sub-basin in order to evaluate and compare the accuracy of the models in predicting several common hydrological and chemical parameters. Batchawana sub-basin is the headwater subbasin of the Turkey Lakes Watershed and is approximately 2.06 km 2 in surface area and receives approximately 1212 mm of precipitation annually. A detailed description of the Turkey Lakes Watershed and its subbasins has been given in Jeffries and Semkin (1982). The ILWAS, ETD and TMWAM models operate on a daily time step but the RAINS model operates on a monthly time step. Also, as not all of the models are capable of simulating all of the major cations and
anions, pH, alkalinity, and sulfate, which are common to all the models, are discussed for comparison purposes• All of the models were applied for the period 1 Jan. 1981-31 Dec. 1984 for data measured at station SI, shown in Fig. 1. The models were initially calibrated for the period 1 Jan. 1981-31 Dec. 1982. They were then run for the final three years of data without further calibration• The models were run using the same basic physical and chemical input constants and parameter values as well as meteorological input data.. Every effort was made to calibrate the different models to the best degree of accuracy, and, where available, the model developer's procedure for model calibration was followed. The models were all calibrated by trial and error rather than by an automatic numerical optimization procedure. The model parameters were adjusted to minimize the errors between the observed and model calculated outputs, based on the principle of least squares. Additional statistics are used to evaluate model performance. Linear and rank correlation coefficients and regression slopes based on a statistical comparison
248
W. G. Booty, A. G. Bobba, D. C. L. Lain, D. S. Jeffries
of computed and observed values are also used to indicate model performance. In the rank correlation, the correlation coefficient (r) and the regression slope (s) are derived from the ordered ranks of both the measured and computed data and are good indicators of episodic coincidences. A correlation coefficient and slope close to 1.0 are indicative of good model performance (Bobba & Lam, 1988). The coefficient of efficiency (E) (Bobba et al., 1988) is a measure of the improvement on the variance not explainable by the mean value, but explainable by the particular model which produces the computed values (qc). The coefficient quantifies the relative improvement of the predictions over a benchmark model (i.e. the mean observed value). The coefficient of efficiency (E) (Nash & Sutcliffe, 1970; Bobba & Lam, 1988) is defined as:
E~
predicting monthly outflows from the watershed. The hydrologic submodel used in the ETD model is intermediate in complexity between the ILWAS and TMWAM and RAINS models. However, for this particular application, it shows the largest deviation from the measured values for both the monthly and cumulative flows, as reflected by all four model performance statistics. Considering each model performance statistic separately, it can be seen that, with respect to the linear correlation coefficient, the RAINS and TMWAM models performed the best with values of r -- 0-987 and 0.986, respectively, followed closely by the ILWAS model with r -- 0.967, and the ETD model was the least precise with r = 0.896. The rank correlation coefficient indicates how well the models simulate the timing and magnitude of the peaks of flows as compared to the linear correlation coefficient. In this case, the TMWAM, ILWAS, and RAINS models perform equally well, with R = 0-988, 0.985, and 0.984, respectively. Considering the coefficient of efficiency, the RAINS and TMWAM models perform the best with respect to the mean-asmodel, with E = 0.970 and 0.961, respectively, followed closely by the ILWAS model with R = 0.915, and finally the ETD model with R -- 0.762. The models show a significant range in their abilities to simulate the monthly and seasonal changes represented by the observed data for outlet chemistry. The model predicted versus observed outlet stream pH, alkalinity, and sulfate data are presented in Figs 3-5, respectively, and the associated statistical data on model performance with respect to these chemical parameters are shown in Table 2. Visually, it can be seen in Figs 3 and 4 that the RAINS model tends to overpredict the fluctuations of observed pH and alkalinity. The ETD model tends to underpredict the fluctuations of observed
[~(qo,--~lo)2--~(qo,-2]-q¢,) i=1
i=1
E (qo,-
(1) 2
i=1
where qo is the mean of the observed values. Any positive value of E can be interpreted as improvement over the mean-as-model, the closer to + 1 the better. The model simulated versus observed monthly watershed outflow data are presented in Fig. 2 and the associated model performance satatistical data are presented in Table 2. The TMWAM and RAINS models show very similar simulated flows and have almost identical statistical performance values. The ILWAS model, which has the most sophisticated hydrologic submodel, performed equally as well as the TMWAM and RAINS models, except with respect to its coefficient of efficiency, in 600
500
400 ~1<
"
300
,T
200 ,-.-.
,oo 0
///
),.. ........I, -
-
1981 •
Obs.Data
1982 TMV~M
1983 ----
ETD
: ......
1984 RAINS
ILWAS
Fig. 2. Comparison of model simulated versus observed monthly outflow data for Batchawana Lake Watershed at station S1 for the period 1 Jan. 1981-31 Dec. 1984.
Application of watershed acidification models
249
Table 2. Statistical evaluation of model performance for (a) the calibration period (1 Jan. 1981-31 Dec. 1982); (b) period 1 Jan. 1983-31 Dec. 1984; and (c) period 1 Jan. 1981-31 Dec. 1984 Model
Linear correlation coefficient (r)
Rank correlation coefficient (R)
ILWAS TMWAM ETD RAINS
0.967 0.986 0-896 0.987
0.985 0.988 0.892 0.984
ILWAS TMWAM ETD RAINS
0-293 0.438 0.549 0.473
0.613 0-713 0-477 0-177
ILWAS TMWAM ETD RAINS
0-442 0.598 0.781 0.356
0.511 0.533 0.687 0.407
ILWAS TMWAM ETD RAINS
0.488 0.682 0.274 0.126
0.610 0.680 0-233 0.153
Slope (s)
Coefficient of efficiency (E)
Outflow 1.100 0.912 1-181 1.065
0.915 0.961 0.762 0.970
-6-772 0.648 1.616 0.142
-2.965 0.264 0-240 -4.004
0.441 0.851 0.959 0.281
0-275 0.340 0.608 -1.048
0.334 0.819 0.444 -0.172
-4-314 0.423 -0.904 -6-479
pH
A/k
SO,
Values of r and R between 0.8 and 1.0 are considered a strong correlation. A value of s between 0.9 and 1-1 is considered to be a strong correlation. The coefficient of efficiency (E) is a measure of the improvement on the variance. Any positive value of E can be interpreted as improvement over the meanas-model, the closer to + 1 the better. pH but does the best job of simulating the observed alkalinities. The T M W A M model appears to do the best job of simulating the observed pH data and is second to the E T D model in simulating the measured alkalinity data. As seen in Table 2, the T M W A M model sulfate predictions have the highest correlation with the ob-
served data (r = 0-682, R = 0.680, s = 0.819). It is also the only model that has a positive coefficient of efficiency ( E = 0.423). In this case, the other three models' predictions of sulfate concentration in the outlet stream cannot be considered as improvement over the mean-asmodel. In Fig. 5 the observed versus model-predicted data are presented for the concentrations of sulfate at
6.5
6
-r" 5.5
5 ¸
4.5
1981
Oba.Data
1982
- -
TMWAM
1983
---
ETD
....... R A I N S
1984
- -
ILWAS
Fig. 3. Comparison of model simulated versus observed monthly averaged outflow pH values for Batchawana Lake Watershed at station SI for the period 1 Jan. 1981-31 Dec. 1984.
W. G. Booty, A. G. Bobba, D. C. L. Lain, D. S. Jeffries
250 150
,j
/i
100
g~
co
//i
50
<
0
I 1981
Obe.Data
Fig. 4.
I
1982
TMWAM
1983
-
--
ETD
1984
RAINS
- -
ILWAS
Comparison of model simulated versus observed monthly averaged outflow alkalinities for Batchawana Lake Watershed at station S1 for the period 1 Jan. 1981-31 Dec. 1984.
station S I. TMWAM predictions most closely simulate both the magnitude and the timing of the observed sulfate stream concentration peaks, as confirmed by the rank correlation coefficient of 0-680 (Table 2). The ILWAS model also performs well in predicting the timing of the peaks in stream sulfate concentration, with R--0.610, as shown in Table 2. However, the magnitudes of the peaks, especially in the spring of 1983 and 1984, are significantly different from the observed values. As was the case for pH, the ETD model predictions for SO4 show the smallest fluctuations of magnitude over the period of study. However, the timing of the peaks do not coincide very well with those of the observed data, as seen by the low rank correlation coefficient (R = 0-233).
The RAINS model underpredicts the stream sulfate concentrations over the entire period of study. It also fails to predict the timing of the peaks, as verified by the very low rank correlation coefficient shown in Table 2 (R = 0.153).
DISCUSSION The models use different definitions for ANC, as listed in Table 3. This may result in the values being different from those obtained using the strict Gran titration definition for the proton reference level. However, for the pH range (5.5-6.7) observed in the lake outflows, the
200
150
,j
2 100 (o c~ n -I 09
5O
I
1981 Oba.Dat8
Fig. 5.
I
1982 - -
TMWAM
1983 ---
ETD
........ R A I N S
1984 - -
ILWAS
Comparison of model simulated versus observed monthly averaged outflow sulfate concentrations for Batchawana Lake Watershed at station S1 for the period 1 Jan. 1981-31 Dec. 1984.
Application of watershed acidification models Table 3. Model definitions of acid neutralizing capacity (ANC) (Modified from Jenne et al., 1989)
Model
ANC definition
ILWAS
ANC -- [HCO~] + 2[CO 2-] + [OH-] + [H2R'"-] + 2[HR '''2-] + 3[R''a-] + JR'-] + [AIOH 2+] + 2[AI(OH)~] + 3[AI(OH)3°] + 4[AI(OH)g] + 3[AIR''~] + lAIR '2+] + 2[AI(R')~] + 3[AI(R')~ - [H +]
ETD
ANC = [HCO3] + 2[CO~ ] + [OH-] - [H +] + [R' ]
251
structure are responsible for the different levels of performance shown by the models. Many models are calibrated and verified over similar short periods of observation and then are used to make future predictions for longer periods of time (i.e. 10, 20, 50 years). It is clearly advantageous to evaluate the relative uncertainties in model predictions using multiple model predictions for measured data before making long-term forecasts of model responses to changes in such inputs as acid loadings.
TMWAM ANC = ([H ÷] + [Na +] + 2[Ca2+] + 2[Mg2+]) - ([HCO;] + 2[SO;])
RAINS
ANC = [HCO;] - [H +] - [AP÷]
Square brackets indicate molar units of concentration, and R', R", and R'" represent monoprotic, diprotic, and triprotic organic acids, respectively.
differences in the ANC values calculated by ILWAS and T M W A M should be comparable to the analytically defined A N C values. The ANC values presented for the measured data were determined using a Metrohm model 636 Titroprocessor (Jeffries et al., 1988) which detects the titration inflection point from the second derivative of the curve. Rigorous checking of this methodology has shown that the results obtained are equivalent to those determined using the FI Gran function. No correction factors have been applied to any of the model-predicted ANC values used in the comparisons presented in this paper. As it has been shown in numerous studies (Johnson et al., 1969; Jeffries & Semkin, 1982; Bobba et al., 1986; Rustad et al., 1986), watershed stream chemistry is strongly correlated with stream hydrology. It might be assumed that, if the hydrological submodel of the ETD model was improved so that its performance was more similar to those shown in Table 2 for the other three models, its ability to predict stream chemistry parameters such as pH, alkalinity, and sulfate may also be improved. The results of these comparisons of model performance for Batchawana Lake Watershed outflow are presented only as an example of how well the models perform for this one site, at one point, and based on monthly averaged data. If we were to apply the models to another period of observations, and/or to different watershed parameters, the relative performance of the models may be expected to be different than those presented in this study. However, the main point is to show that basic watershed parameters, such as basin outflow, outflow pH, alkalinity, and sulfate are simulated to different degrees of accuracy by the four different models considered in this study. A great deal more effort is still required in order to fully test and validate the models, and more importantly, to help to understand which key processes or differences in model
ACKNOWLEDGMENTS The assistance of R. Semkin in the acquisition of watershed data is gratefully acknowledged. The authors are also very grateful to H. M. Seip, N. Christophersen, S. A. Gherini, R. K. Munson, M. M. Lang, J. L. Schnoor, N. P. Nickolaidas and J. Kamari for providing the various model codes and for useful discussions.
REFERENCES
Arp, P. A. & Ramnarine, S. (1983). Verifying the effects of acid precipitation on soil leachates--A comparison between published records and model predictions. Ecol. Modelling, 19, 119-38. Beck, M. B., Drummond, D., Kleissen, F. M., Langan, S. J., Wheater, H. S. & Whitehead, P. G. (1987). Surface water acidification: The Birkenes data revisited. In Systems Analysis in Water Quality Management, ed. M. B. Beck, IAWPRC/Pergamon, New York, pp. 133-50. Bobba, A. G. & Lam, D. C. L. (1988). Application of a hydrological model to the Turkey Lakes watershed. Can. J. Fish. Aquat. Sci., 45, 81-7. Bobba, A. G. & Lam, D. C. L. (1989). Application of a hydrological model to acidified watersheds: A study on Mersey River and Moosepit Brook, Nova Scotia. Water, Air & Soil Poll., 46, 261-75. Bobba, A. G. & Lam, D. C. L. (1990). Hydrological modelling of acidified Canadian watersheds. Ecol. Modelling, 50, 532. Bobba, A. G., Lam, D. C. L. & De Grosbois, E. (1986). Application of hydrological model to acidified watersheds: A study on Harp Lake catchments. Hydrol. Sci. Technol. Short Pap., 2, 37-44. Booty, W. G. (1983). Watershed acidification model and the soil acid neutralization capacity concept, PhD dissertation, McMaster University, Hamilton, Ontario. Booty, W. G. & Kramer, J. R. (1984). Sensitivity analysis of a watershed acidification model. Philos. Trans. R. Soc. London, B305, 441-9. Chen, C. W., Gherini, S. A., Hodson, R. J. M. & Dean, J. D. (1983). The ILWAS study: Model principles and applications procedures. EPRI EA-3221, Project 1109-5, Vol. 1, Electric Power Research Institute, Palo Alto, CA. Christophersen, N., Seip, H. M. & Wright, R. F., (1982). A model for streamwater chemistry at Birkenes, Norway. Water Resour. Res., 18, 977-96. Cosby, B. J., Hornberger, G. M., Galloway, J. N. & Wright, R. F. (1985a). Modelling the effects of acid deposition: Assessment of a lumped parameter model of soil water and stream chemistry. Water Resour. Res., 21, 51-63.
252
W. G. Booty, A. G. Bobba, D. C. L. Lam, D. S. Jeffries
Cosby, B. J., Hornberger, G. M., Galloway, J. N. & Wright, R. F. (1985b). Time scales of catchment acidification: A quantitative model for estimating freshwater acidification. Environ. Sct Technol., 19, 1144-9. Davis, J. E. & Goldstein, R. A. (1988). Simulated response of an acidic Adirondack lake-watershed to various liming mitigation strategies. Water Resour. Res., 24, 525-32. Jeffries, D. S. & Semkin, R. G. (1982). Basin description and information pertinent to mass balance studies of the Turkey Lakes watershed. Turkey Lakes Watershed Rep. TLW-8201, Can. Centre Inland Waters, Burlington, Ontario. Jeffries, D. S. & Semkin, R. G. (1983). Changes in snowpack, stream, and lake chemistry during snowmelt in the Turkey Lakes Watershed. VDI-Berichte, 500, 377-86. Jeffries, D. S., Semkin, R. G., Neureuther, R. & Seymour, M. (1988). Ion mass budgets for lakes in the Turkey Lakes watershed, June 1981-May, 1983. Can. J. Fish. Aquat. Sci., 45, 47-58. Jenne, E. A. et al. (1989). An evaluation and analysis of three dynamic watershed acidification codes (MAGIC, ETD, and ILWAS). US Environmental Protection Agency, Washington, DC, PNL-6687, UC-11. Johnson, N. M., Likens, G. E., Bormann, F. H., Fisher, D. W. & Pierce, R. S. (1969). A working model for the variation of stream water chemistry at the Hubbard Brook Experimental Forest, New Hampshire. Water Resour. Res., 5, 1353-63. Kamari, J. (1985). A model for analyzing lake water acidification on a regional scale. Report no. CP-85-48. International Institute for Applied Systems Analysis, 2361 Laxenberg, Austria.
Lam, D. C. L., Bobba, A. G., Jeffries, D. S. & Craig, D. (1988). Modelling stream chemistry for the Turkey Lakes Watershed: Comparison with 1981-84 data. Can. J. Fish. Aquat. Sci., 45, 72-80. Lam, D. C. L., Bobba, A. G., Bourbonniere, R. A., Howell, G. D. & Thompson, M. E. (1989). Modelling organic and inorganic acidity in two Nova Scotia rivers. Water, Air & Soil Poll., 46, 277-87. Nash, J. E. & Sutcliffe, J. V. (1970). Riverflow forecasting through conceptual models. 1. A discussion of principles. J. Hydrol., 10, 282-90. Nikolaidas, N. P., Rajaram, H., Schnoor, J. L. & Geogakakos, K. P. (1988). A generalized soft water acidification model. Water Resour. Res., 24, 1983-96. Ormerod, S. J., Weatherly, N. S. & Whitehead, P. G. (1988). Temporal patterns of ecological change during the acidification and recovery of streams throughout Wales according to a hydrochemical model. In Regional Acidification Models, ed. J. Kamarai. Springer-Verlag, New York, pp. 167-84. RMCC (Research Monitoring Coordinating Committee) (1990). Canadian Long-Range Transport of Air Pollutants and Acid Deposition Assessment Report. Volumes !-8, Environment Canada, Ottawa. Rustad, S., Christophersen, N., Seip, H. M. & Dillon, P. J. (1986). Model for streamwater chemistry of a tributary to Harp Lake, Ontario. Can. J. Fish. Aquat. Sci., 43, 625-33. Schnoor, J. L., Palmer, W. D., Jr & Glass, G. E. (1984). Modeling impacts of acid precipitation for northeastern Minnesota. In Modeling of Total Acid Precipitation Impact, ed. J. L. Schnoor. Butterworth, Boston, MA, pp. 155-73.